Dietary intake, Nutritional status, and Health outcomes among Vegan, Vegetarian and Omnivore families: results from the observational study
Statistical report - diet prediction (elastic net)
Authors and affiliations
Marina Heniková1,2, Anna Ouřadová1, Eliška Selinger1,3, Filip Tichanek4, Petra Polakovičová4, Dana Hrnčířová2, Pavel Dlouhý2, Martin Světnička5, Eva El-Lababidi5, Jana Potočková1, Tilman Kühn6, Monika Cahová4, Jan Gojda1
1 Department of Internal Medicine, Kralovske Vinohrady University Hospital and Third Faculty of Medicine, Charles University, Prague, Czech Republic.
2 Department of Hygiene, Third Faculty of Medicine, Charles University, Prague, Czech Republic.
3 National Health Institute, Prague, Czech Republic.
4 Institute for Clinical and Experimental Medicine, Prague, Czech Republic.
5 Department of Pediatrics, Kralovske Vinohrady University Hospital and Third Faculty of Medicine, Charles University, Prague, Czech Republic.
6 Department of Epidemiology, MedUni, Vienna, Austria.
This is a statistical report of the study currenlty under review in the Communications Medicine journal.
When using this code or data, cite the original publication:
TO BE ADDED
BibTex citation for the original publication:
TO BE ADDED
Original GitHub repository: https://github.com/filip-tichanek/kompas_clinical
Statistical reports can be found on the reports hub.
Data analysis is described in detail in the statistical methods report.
1 Introduction
This project is designed to evaluate and compare clinical outcomes across three distinct dietary strategy groups:
- Vegans
- Vegetarians
- Omnivores
The dataset includes both adults and children, with data clustered within families.
1.1 Main Questions
The study addresses the following key questions:
Q1. Do clinical outcomes vary significantly across different diet strategies?
Q2. Beyond diet group, which factors (e.g., sex, age, breastfeeding status for children, or supplementation when applicable) most strongly influence clinical outcomes? How correlated (“clustered”) are these characteristics within the same family?
Q3. Could the clinical characteristics effectively discriminate between different diet groups?
1.2 Statistical Methods
For full methodological details, see this report. In brief:
Robust linear mixed-effects models (rLME) were used to estimate adjusted differences between diet groups (Q1) and assess the importance of other variables (Q2), including how much clinical characteristics tend to cluster within families. Covariates included age, sex, breastfeeding status for children, and relevant supplementation factors where applicable.
Elastic net logistic regression was employed to answer Q3, evaluating whether clinical characteristics provide a strong overall signal distinguishing between diet groups, incorporating a predictive perspective.
All analyses were conducted separately for adults and children.
Open code
# Import initiation file
source('r/368_initiation.R')2 Data wrangling
2.1 All children
Open code
skimr::skim(dat_child_all)| Name | dat_child_all |
| Number of rows | 142 |
| Number of columns | 75 |
| _______________________ | |
| Column type frequency: | |
| factor | 4 |
| numeric | 71 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| ID | 0 | 1 | FALSE | 142 | K10: 1, K10: 1, K10: 1, K11: 1 |
| FAM | 0 | 1 | FALSE | 95 | R25: 3, R41: 3, R55: 3, R70: 3 |
| SEX | 0 | 1 | FALSE | 2 | F: 73, M: 69 |
| GRP | 0 | 1 | FALSE | 3 | VN: 61, OM: 44, VG: 37 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| aAGE | 0 | 1.00 | 3.34 | 2.46 | 0.18 | 1.34 | 2.74 | 4.88 | 12.04 | ▇▅▃▁▁ |
| aBreastFeed_full | 0 | 1.00 | 6.04 | 2.10 | 0.00 | 6.00 | 6.00 | 6.00 | 18.00 | ▁▇▁▁▁ |
| aBreastFeed_total | 13 | 0.91 | 15.90 | 10.14 | 0.00 | 9.21 | 14.00 | 19.00 | 60.00 | ▇▇▂▁▁ |
| aBreastFeed_full_stopped | 0 | 1.00 | 0.92 | 0.27 | 0.00 | 1.00 | 1.00 | 1.00 | 1.00 | ▁▁▁▁▇ |
| aBreastFeed_total_stopped | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| dev_delay | 0 | 1.00 | 0.02 | 0.14 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aMASS_Perc | 0 | 1.00 | 46.85 | 28.80 | 0.00 | 21.00 | 49.00 | 68.50 | 99.00 | ▇▆▃▇▅ |
| aHEIGHT_Perc | 0 | 1.00 | 45.35 | 27.47 | 0.00 | 19.50 | 44.00 | 68.75 | 100.00 | ▇▇▅▇▃ |
| aBMI_PERC | 0 | 1.00 | 50.50 | 24.96 | 0.00 | 34.25 | 50.00 | 68.00 | 100.00 | ▃▅▇▅▃ |
| aM_per_H_PERC | 0 | 1.00 | 49.98 | 24.85 | 0.00 | 33.00 | 51.00 | 66.75 | 99.00 | ▃▃▇▆▃ |
| aGLY | 10 | 0.93 | 4.45 | 0.59 | 3.26 | 4.16 | 4.42 | 4.72 | 7.49 | ▂▇▁▁▁ |
| aTC | 10 | 0.93 | 3.78 | 0.66 | 2.21 | 3.30 | 3.87 | 4.21 | 5.56 | ▂▆▇▅▁ |
| aHDL | 10 | 0.93 | 1.26 | 0.33 | 0.66 | 1.01 | 1.25 | 1.42 | 2.16 | ▆▇▇▂▂ |
| aLDL | 10 | 0.93 | 2.05 | 0.57 | 0.24 | 1.69 | 2.02 | 2.41 | 3.48 | ▁▃▇▆▂ |
| aTG | 10 | 0.93 | 1.06 | 0.60 | 0.30 | 0.66 | 0.92 | 1.22 | 3.57 | ▇▅▁▁▁ |
| aCa | 9 | 0.94 | 2.57 | 0.12 | 2.35 | 2.48 | 2.56 | 2.65 | 2.88 | ▃▇▆▅▁ |
| aP | 9 | 0.94 | 1.66 | 0.18 | 1.30 | 1.54 | 1.67 | 1.79 | 2.10 | ▃▆▇▅▁ |
| aMg | 9 | 0.94 | 0.86 | 0.07 | 0.71 | 0.81 | 0.85 | 0.90 | 1.17 | ▃▇▃▁▁ |
| aSe | 11 | 0.92 | 0.79 | 0.28 | 0.27 | 0.60 | 0.77 | 0.94 | 2.02 | ▅▇▃▁▁ |
| aZn | 11 | 0.92 | 11.56 | 2.41 | 7.00 | 10.15 | 11.30 | 12.80 | 21.30 | ▃▇▃▁▁ |
| aFE | 9 | 0.94 | 13.37 | 6.94 | 1.60 | 8.30 | 12.60 | 17.80 | 32.70 | ▆▇▇▂▁ |
| aVKFE | 10 | 0.93 | 71.66 | 9.12 | 48.00 | 66.00 | 72.00 | 77.00 | 103.00 | ▂▆▇▂▁ |
| aFERR | 9 | 0.94 | 18.62 | 15.03 | 3.10 | 10.40 | 15.20 | 21.70 | 133.70 | ▇▁▁▁▁ |
| aTRF | 9 | 0.94 | 2.84 | 0.36 | 1.90 | 2.60 | 2.85 | 3.05 | 4.10 | ▂▆▇▂▁ |
| aSATTRF | 9 | 0.94 | 18.90 | 9.89 | 2.10 | 11.60 | 18.10 | 24.90 | 43.70 | ▅▇▆▂▂ |
| aTRFINDEX | 9 | 0.94 | 1.55 | 0.74 | 0.59 | 1.13 | 1.38 | 1.70 | 5.41 | ▇▃▁▁▁ |
| aSTRF | 9 | 0.94 | 1.67 | 0.33 | 1.08 | 1.43 | 1.62 | 1.81 | 2.98 | ▅▇▂▁▁ |
| aHGB | 14 | 0.90 | 121.94 | 9.48 | 99.00 | 115.00 | 121.00 | 127.25 | 157.00 | ▁▇▆▂▁ |
| aMCV | 14 | 0.90 | 78.81 | 3.83 | 68.20 | 76.77 | 79.45 | 81.10 | 88.30 | ▂▃▇▇▁ |
| aPTH | 10 | 0.93 | 2.88 | 1.25 | 0.70 | 1.90 | 2.80 | 3.60 | 7.80 | ▆▇▃▁▁ |
| aCros | 11 | 0.92 | 1.26 | 0.31 | 0.57 | 1.08 | 1.23 | 1.42 | 2.71 | ▃▇▃▁▁ |
| aP1NP | 11 | 0.92 | 720.87 | 294.64 | 215.00 | 482.00 | 631.00 | 959.50 | 1200.99 | ▂▇▅▂▅ |
| aUI | 30 | 0.79 | 176.10 | 142.86 | 6.00 | 86.95 | 138.00 | 228.50 | 911.00 | ▇▃▁▁▁ |
| aUREA | 9 | 0.94 | 4.21 | 1.24 | 1.80 | 3.40 | 4.20 | 5.10 | 7.70 | ▅▇▇▃▁ |
| aCREA | 9 | 0.94 | 25.99 | 7.50 | 15.00 | 20.00 | 25.00 | 31.00 | 51.00 | ▇▅▅▁▁ |
| aUA | 9 | 0.94 | 226.97 | 48.32 | 78.00 | 197.00 | 224.00 | 257.00 | 391.00 | ▁▅▇▃▁ |
| aVIT_AKTB12 | 11 | 0.92 | 127.83 | 82.79 | 31.50 | 77.75 | 108.70 | 144.05 | 493.90 | ▇▃▁▁▁ |
| aHCY | 23 | 0.84 | 8.69 | 3.09 | 3.90 | 6.85 | 8.10 | 9.75 | 25.80 | ▇▆▁▁▁ |
| aMMA | 19 | 0.87 | 286.37 | 354.72 | 67.30 | 148.95 | 187.70 | 258.25 | 2980.40 | ▇▁▁▁▁ |
| aVIT_D | 9 | 0.94 | 92.73 | 30.42 | 36.90 | 69.80 | 87.60 | 110.20 | 194.20 | ▅▇▅▂▁ |
| aFOLAT | 10 | 0.93 | 16.41 | 4.66 | 4.50 | 13.05 | 16.90 | 19.65 | 24.04 | ▁▅▇▇▆ |
| aIGF1 | 20 | 0.86 | 110.37 | 63.87 | 24.60 | 63.60 | 92.35 | 136.98 | 318.00 | ▇▆▃▁▁ |
| aUr_Ca | 18 | 0.87 | 1.55 | 1.55 | 0.25 | 0.37 | 1.02 | 2.26 | 7.41 | ▇▂▁▁▁ |
| aCa_per_Krea | 18 | 0.87 | 0.48 | 0.51 | 0.03 | 0.12 | 0.31 | 0.64 | 2.75 | ▇▂▁▁▁ |
| aI_per_Krea | 31 | 0.78 | 553.90 | 748.54 | 60.27 | 189.31 | 351.80 | 563.53 | 5887.55 | ▇▁▁▁▁ |
| aP_per_Krea | 18 | 0.87 | 3.54 | 2.13 | 0.14 | 1.71 | 3.14 | 4.73 | 9.24 | ▇▇▆▂▂ |
| aAL_child | 0 | 1.00 | 0.12 | 0.33 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aBREAKS | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_VEG1 | 0 | 1.00 | 0.05 | 0.22 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSup_B12 | 0 | 1.00 | 0.50 | 0.50 | 0.00 | 0.00 | 0.50 | 1.00 | 1.00 | ▇▁▁▁▇ |
| aSUP_FOL | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_vitA | 0 | 1.00 | 0.01 | 0.12 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_vitB1 | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_vit.B5 | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_D | 0 | 1.00 | 0.72 | 0.45 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| aSUP_Mg | 0 | 1.00 | 0.05 | 0.22 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Zn | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Se | 0 | 1.00 | 0.02 | 0.14 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Ca | 0 | 1.00 | 0.04 | 0.18 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Fe | 0 | 1.00 | 0.07 | 0.26 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Iod | 0 | 1.00 | 0.06 | 0.24 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_Ѡ3 | 0 | 1.00 | 0.40 | 0.49 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▆ |
| aSUP_CHLO | 0 | 1.00 | 0.02 | 0.14 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_ALG | 0 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | ▁▁▇▁▁ |
| aSUP_GB | 0 | 1.00 | 0.01 | 0.08 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_FORT | 0 | 1.00 | 0.01 | 0.08 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_PROB | 0 | 1.00 | 0.06 | 0.24 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| aSUP_OTH | 0 | 1.00 | 0.21 | 0.41 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| aBiW | 2 | 0.99 | 3336.26 | 541.50 | 2000.00 | 3047.50 | 3400.00 | 3657.50 | 5000.00 | ▂▅▇▃▁ |
| log2_age | 0 | 1.00 | 1.29 | 1.25 | -2.49 | 0.43 | 1.45 | 2.29 | 3.59 | ▁▃▆▇▅ |
| aBreastFeed_full_duration | 0 | 1.00 | 5.64 | 2.61 | 0.00 | 5.50 | 6.00 | 6.00 | 18.00 | ▁▇▁▁▁ |
2.1.1 Wrangle data
Open code
for (i in 1:ncol(dat_child_all)) {
names(dat_child_all)[i] <- str_replace_all(
names(dat_child_all)[i], c("^a" = "")
)
}
dat_child_all <- dat_child_all %>%
rename(
FOLATE = FOLAT,
TIBC = VKFE,
HOLOTC = VIT_AKTB12,
CTX = Cros,
uCCR = Ca_per_Krea,
uICR = I_per_Krea,
uPCR = P_per_Krea,
uCa = Ur_Ca
) %>%
mutate(
log2_TG = log2(TG),
log2_Se = log2(Se),
log2_Zn = log2(Zn),
log2_FERR = log2(FERR),
log2_TRFINDEX = log2(TRFINDEX),
log2_STRF = log2(STRF),
log2_CTX = log2(CTX),
log2_UI = log2(UI),
log2_CREA = log2(CREA),
log2_HOLOTC = log2(HOLOTC),
log2_HCY = log2(HCY),
log2_MMA = log2(MMA),
log2_VIT_D = log2(VIT_D),
log2_IGF1 = log2(IGF1),
log2_uCa = log2(uCa),
log2_uCCR = log2(uCCR),
log2_uICR = log2(uICR)
) %>%
mutate(
male = if_else(SEX == "M", 1, 0),
GRP = factor(GRP),
FAM = factor(FAM)
) %>%
select(
## remove irrelevant variables not included to any inference and prediction
-c(SEX, BreastFeed_full, BreastFeed_total, ID, dev_delay, BREAKS, AGE),
## remove information about supplementation
-c(SUP_VEG1:SUP_OTH),
## remove non-log2-transformed version when log2(variable) is included
-c(TG, Se, Zn, FERR, TRFINDEX,
STRF, CTX, UI, CREA, HOLOTC, HCY,
MMA, VIT_D, IGF1, uCa, uCCR, uICR
)
) %>%
select(FAM, GRP, log2_age, male, BiW,
BreastFeed_full_duration, BreastFeed_full_stopped, BreastFeed_total_stopped,
everything()
)2.1.2 Missing values and standardization
Open code
skimr::skim(dat_child_all)| Name | dat_child_all |
| Number of rows | 142 |
| Number of columns | 49 |
| _______________________ | |
| Column type frequency: | |
| factor | 2 |
| numeric | 47 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| FAM | 0 | 1 | FALSE | 95 | R25: 3, R41: 3, R55: 3, R70: 3 |
| GRP | 0 | 1 | FALSE | 3 | VN: 61, OM: 44, VG: 37 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| log2_age | 0 | 1.00 | 1.29 | 1.25 | -2.49 | 0.43 | 1.45 | 2.29 | 3.59 | ▁▃▆▇▅ |
| male | 0 | 1.00 | 0.49 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| BiW | 2 | 0.99 | 3336.26 | 541.50 | 2000.00 | 3047.50 | 3400.00 | 3657.50 | 5000.00 | ▂▅▇▃▁ |
| BreastFeed_full_duration | 0 | 1.00 | 5.64 | 2.61 | 0.00 | 5.50 | 6.00 | 6.00 | 18.00 | ▁▇▁▁▁ |
| BreastFeed_full_stopped | 0 | 1.00 | 0.92 | 0.27 | 0.00 | 1.00 | 1.00 | 1.00 | 1.00 | ▁▁▁▁▇ |
| BreastFeed_total_stopped | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| MASS_Perc | 0 | 1.00 | 46.85 | 28.80 | 0.00 | 21.00 | 49.00 | 68.50 | 99.00 | ▇▆▃▇▅ |
| HEIGHT_Perc | 0 | 1.00 | 45.35 | 27.47 | 0.00 | 19.50 | 44.00 | 68.75 | 100.00 | ▇▇▅▇▃ |
| BMI_PERC | 0 | 1.00 | 50.50 | 24.96 | 0.00 | 34.25 | 50.00 | 68.00 | 100.00 | ▃▅▇▅▃ |
| M_per_H_PERC | 0 | 1.00 | 49.98 | 24.85 | 0.00 | 33.00 | 51.00 | 66.75 | 99.00 | ▃▃▇▆▃ |
| GLY | 10 | 0.93 | 4.45 | 0.59 | 3.26 | 4.16 | 4.42 | 4.72 | 7.49 | ▂▇▁▁▁ |
| TC | 10 | 0.93 | 3.78 | 0.66 | 2.21 | 3.30 | 3.87 | 4.21 | 5.56 | ▂▆▇▅▁ |
| HDL | 10 | 0.93 | 1.26 | 0.33 | 0.66 | 1.01 | 1.25 | 1.42 | 2.16 | ▆▇▇▂▂ |
| LDL | 10 | 0.93 | 2.05 | 0.57 | 0.24 | 1.69 | 2.02 | 2.41 | 3.48 | ▁▃▇▆▂ |
| Ca | 9 | 0.94 | 2.57 | 0.12 | 2.35 | 2.48 | 2.56 | 2.65 | 2.88 | ▃▇▆▅▁ |
| P | 9 | 0.94 | 1.66 | 0.18 | 1.30 | 1.54 | 1.67 | 1.79 | 2.10 | ▃▆▇▅▁ |
| Mg | 9 | 0.94 | 0.86 | 0.07 | 0.71 | 0.81 | 0.85 | 0.90 | 1.17 | ▃▇▃▁▁ |
| FE | 9 | 0.94 | 13.37 | 6.94 | 1.60 | 8.30 | 12.60 | 17.80 | 32.70 | ▆▇▇▂▁ |
| TIBC | 10 | 0.93 | 71.66 | 9.12 | 48.00 | 66.00 | 72.00 | 77.00 | 103.00 | ▂▆▇▂▁ |
| TRF | 9 | 0.94 | 2.84 | 0.36 | 1.90 | 2.60 | 2.85 | 3.05 | 4.10 | ▂▆▇▂▁ |
| SATTRF | 9 | 0.94 | 18.90 | 9.89 | 2.10 | 11.60 | 18.10 | 24.90 | 43.70 | ▅▇▆▂▂ |
| HGB | 14 | 0.90 | 121.94 | 9.48 | 99.00 | 115.00 | 121.00 | 127.25 | 157.00 | ▁▇▆▂▁ |
| MCV | 14 | 0.90 | 78.81 | 3.83 | 68.20 | 76.77 | 79.45 | 81.10 | 88.30 | ▂▃▇▇▁ |
| PTH | 10 | 0.93 | 2.88 | 1.25 | 0.70 | 1.90 | 2.80 | 3.60 | 7.80 | ▆▇▃▁▁ |
| P1NP | 11 | 0.92 | 720.87 | 294.64 | 215.00 | 482.00 | 631.00 | 959.50 | 1200.99 | ▂▇▅▂▅ |
| UREA | 9 | 0.94 | 4.21 | 1.24 | 1.80 | 3.40 | 4.20 | 5.10 | 7.70 | ▅▇▇▃▁ |
| UA | 9 | 0.94 | 226.97 | 48.32 | 78.00 | 197.00 | 224.00 | 257.00 | 391.00 | ▁▅▇▃▁ |
| FOLATE | 10 | 0.93 | 16.41 | 4.66 | 4.50 | 13.05 | 16.90 | 19.65 | 24.04 | ▁▅▇▇▆ |
| uPCR | 18 | 0.87 | 3.54 | 2.13 | 0.14 | 1.71 | 3.14 | 4.73 | 9.24 | ▇▇▆▂▂ |
| AL_child | 0 | 1.00 | 0.12 | 0.33 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| log2_TG | 10 | 0.93 | -0.10 | 0.71 | -1.74 | -0.60 | -0.12 | 0.29 | 1.84 | ▂▆▇▃▁ |
| log2_Se | 11 | 0.92 | -0.43 | 0.51 | -1.89 | -0.74 | -0.38 | -0.10 | 1.01 | ▁▃▇▅▁ |
| log2_Zn | 11 | 0.92 | 3.50 | 0.29 | 2.81 | 3.34 | 3.50 | 3.68 | 4.41 | ▂▆▇▂▁ |
| log2_FERR | 9 | 0.94 | 3.92 | 0.90 | 1.63 | 3.38 | 3.93 | 4.44 | 7.06 | ▂▆▇▂▁ |
| log2_TRFINDEX | 9 | 0.94 | 0.52 | 0.54 | -0.76 | 0.18 | 0.46 | 0.77 | 2.44 | ▂▇▆▁▁ |
| log2_STRF | 9 | 0.94 | 0.72 | 0.27 | 0.11 | 0.52 | 0.70 | 0.86 | 1.58 | ▂▇▇▂▁ |
| log2_CTX | 11 | 0.92 | 0.29 | 0.36 | -0.81 | 0.11 | 0.30 | 0.50 | 1.44 | ▁▃▇▃▁ |
| log2_UI | 30 | 0.79 | 7.03 | 1.25 | 2.58 | 6.44 | 7.11 | 7.84 | 9.83 | ▁▁▆▇▂ |
| log2_CREA | 9 | 0.94 | 4.64 | 0.40 | 3.91 | 4.32 | 4.64 | 4.95 | 5.67 | ▇▇▇▆▁ |
| log2_HOLOTC | 11 | 0.92 | 6.76 | 0.81 | 4.98 | 6.28 | 6.76 | 7.17 | 8.95 | ▂▆▇▃▁ |
| log2_HCY | 23 | 0.84 | 3.05 | 0.45 | 1.96 | 2.78 | 3.02 | 3.29 | 4.69 | ▁▇▆▁▁ |
| log2_MMA | 19 | 0.87 | 7.75 | 0.91 | 6.07 | 7.22 | 7.55 | 8.01 | 11.54 | ▂▇▁▁▁ |
| log2_VIT_D | 9 | 0.94 | 6.46 | 0.46 | 5.21 | 6.13 | 6.45 | 6.78 | 7.60 | ▁▅▇▅▁ |
| log2_IGF1 | 20 | 0.86 | 6.56 | 0.83 | 4.62 | 5.99 | 6.53 | 7.10 | 8.31 | ▂▆▇▇▃ |
| log2_uCa | 18 | 0.87 | -0.06 | 1.47 | -2.00 | -1.43 | 0.04 | 1.18 | 2.89 | ▇▅▅▆▂ |
| log2_uCCR | 18 | 0.87 | -1.78 | 1.50 | -5.14 | -3.00 | -1.68 | -0.65 | 1.46 | ▂▆▇▆▂ |
| log2_uICR | 31 | 0.78 | 8.46 | 1.28 | 5.91 | 7.56 | 8.46 | 9.14 | 12.52 | ▃▇▇▂▁ |
Removing features with over 25% NA values: all variables are below this threshold, no will be excluded
2.1.2.1 Multiple imputation
Prepare
Open code
dat_child_all_toimp <- dat_child_all %>% select(-FAM)
init <- mice(dat_child_all_toimp, maxit = 0)
init$method
## GRP log2_age male
## "" "" ""
## BiW BreastFeed_full_duration BreastFeed_full_stopped
## "pmm" "" ""
## BreastFeed_total_stopped MASS_Perc HEIGHT_Perc
## "" "" ""
## BMI_PERC M_per_H_PERC GLY
## "" "" "pmm"
## TC HDL LDL
## "pmm" "pmm" "pmm"
## Ca P Mg
## "pmm" "pmm" "pmm"
## FE TIBC TRF
## "pmm" "pmm" "pmm"
## SATTRF HGB MCV
## "pmm" "pmm" "pmm"
## PTH P1NP UREA
## "pmm" "pmm" "pmm"
## UA FOLATE uPCR
## "pmm" "pmm" "pmm"
## AL_child log2_TG log2_Se
## "" "pmm" "pmm"
## log2_Zn log2_FERR log2_TRFINDEX
## "pmm" "pmm" "pmm"
## log2_STRF log2_CTX log2_UI
## "pmm" "pmm" "pmm"
## log2_CREA log2_HOLOTC log2_HCY
## "pmm" "pmm" "pmm"
## log2_MMA log2_VIT_D log2_IGF1
## "pmm" "pmm" "pmm"
## log2_uCa log2_uCCR log2_uICR
## "pmm" "pmm" "pmm"
data_imputed <- run(
expr = mice(
dat_child_all_toimp,
method = init$method,
m = 1
),
path = "gitignore/run/data_child_all_imputed"
)
dat_child_all_imputed <- data.frame(
FAM = dat_child_all$FAM,
complete(data_imputed, 'all')[1]$`1`
)2.1.2.2 Rescaling
Open code
dat_child_all_imputed <- dat_child_all_imputed %>%
mutate(
across(
where(~ is.numeric(.) && length(unique(.)) > 2),
arm::rescale
)
)
summary(dat_child_all_imputed)
## FAM GRP log2_age male BiW
## R25 : 3 OM:44 Min. :-1.51597 Min. :0.0000 Min. :-1.24365
## R41 : 3 VG:37 1st Qu.:-0.34675 1st Qu.:0.0000 1st Qu.:-0.26830
## R55 : 3 VN:61 Median : 0.06553 Median :0.0000 Median : 0.06022
## R70 : 3 Mean : 0.00000 Mean :0.4859 Mean : 0.00000
## R91 : 3 3rd Qu.: 0.39903 3rd Qu.:1.0000 3rd Qu.: 0.30846
## R11 : 2 Max. : 0.92159 Max. :1.0000 Max. : 1.54041
## (Other):125
## BreastFeed_full_duration BreastFeed_full_stopped BreastFeed_total_stopped
## Min. :-1.08165 Min. :0.0000 Min. :0.0000
## 1st Qu.:-0.02701 1st Qu.:1.0000 1st Qu.:0.0000
## Median : 0.06887 Median :1.0000 Median :1.0000
## Mean : 0.00000 Mean :0.9225 Mean :0.6761
## 3rd Qu.: 0.06887 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. : 2.36990 Max. :1.0000 Max. :1.0000
##
## MASS_Perc HEIGHT_Perc BMI_PERC M_per_H_PERC
## Min. :-0.81336 Min. :-0.82527 Min. :-1.01169 Min. :-1.00579
## 1st Qu.:-0.44880 1st Qu.:-0.47038 1st Qu.:-0.32554 1st Qu.:-0.34169
## Median : 0.03729 Median :-0.02448 Median :-0.01002 Median : 0.02055
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.37581 3rd Qu.: 0.42597 3rd Qu.: 0.35058 3rd Qu.: 0.33751
## Max. : 0.90530 Max. : 0.99471 Max. : 0.99165 Max. : 0.98652
##
## GLY TC HDL LDL
## Min. :-0.96932 Min. :-1.1611 Min. :-0.897027 Min. :-1.5591
## 1st Qu.:-0.24605 1st Qu.:-0.3518 1st Qu.:-0.381160 1st Qu.:-0.3186
## Median :-0.02522 Median : 0.0427 Median :-0.008379 Median :-0.0237
## Mean : 0.00000 Mean : 0.0000 Mean : 0.000000 Mean : 0.0000
## 3rd Qu.: 0.23005 3rd Qu.: 0.3271 3rd Qu.: 0.274030 3rd Qu.: 0.3124
## Max. : 2.45859 Max. : 1.2979 Max. : 1.362247 Max. : 1.2515
##
## Ca P Mg FE
## Min. :-0.95206 Min. :-1.04671 Min. :-1.07088 Min. :-0.85319
## 1st Qu.:-0.40979 1st Qu.:-0.35459 1st Qu.:-0.30075 1st Qu.:-0.36807
## Median :-0.03437 Median : 0.02031 Median :-0.02071 Median :-0.05248
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.34105 3rd Qu.: 0.36637 3rd Qu.: 0.31185 3rd Qu.: 0.33165
## Max. : 1.25874 Max. : 1.26036 Max. : 2.14964 Max. : 1.39025
##
## TIBC TRF SATTRF HGB
## Min. :-1.28220 Min. :-1.29718 Min. :-0.86511 Min. :-1.1969
## 1st Qu.:-0.33960 1st Qu.:-0.38031 1st Qu.:-0.37882 1st Qu.:-0.3516
## Median : 0.02925 Median : 0.03263 Median :-0.04314 Median :-0.0346
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.30246 3rd Qu.: 0.29859 3rd Qu.: 0.29509 3rd Qu.: 0.2824
## Max. : 1.72320 Max. : 1.78238 Max. : 1.25873 Max. : 1.8674
##
## MCV PTH P1NP UREA
## Min. :-1.27228 Min. :-0.86222 Min. :-0.8720 Min. :-0.95114
## 1st Qu.:-0.23179 1st Qu.:-0.37956 1st Qu.:-0.4112 1st Qu.:-0.33321
## Median : 0.06504 Median :-0.03767 Median :-0.1493 Median :-0.01939
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.29939 3rd Qu.: 0.30421 3rd Qu.: 0.5204 3rd Qu.: 0.34505
## Max. : 1.23989 Max. : 1.99353 Max. : 0.7624 Max. : 1.43836
##
## UA FOLATE uPCR AL_child
## Min. :-1.4333036 Min. :-1.31482 Min. :-0.7806 Min. :0.0000
## 1st Qu.:-0.3100253 1st Qu.:-0.36494 1st Qu.:-0.4357 1st Qu.:0.0000
## Median :-0.0000698 Median : 0.05861 Median :-0.1030 Median :0.0000
## Mean : 0.0000000 Mean : 0.00000 Mean : 0.0000 Mean :0.1197
## 3rd Qu.: 0.3074059 3rd Qu.: 0.33510 3rd Qu.: 0.2779 3rd Qu.:0.0000
## Max. : 1.6712097 Max. : 0.81004 Max. : 1.2758 Max. :1.0000
##
## log2_TG log2_Se log2_Zn log2_FERR
## Min. :-1.16637 Min. :-1.44955 Min. :-1.176518 Min. :-1.223479
## 1st Qu.:-0.36102 1st Qu.:-0.30729 1st Qu.:-0.259301 1st Qu.:-0.309078
## Median :-0.03236 Median : 0.04957 Median :-0.009786 Median : 0.002261
## Mean : 0.00000 Mean : 0.00000 Mean : 0.000000 Mean : 0.000000
## 3rd Qu.: 0.27003 3rd Qu.: 0.33111 3rd Qu.: 0.322212 3rd Qu.: 0.291327
## Max. : 1.31588 Max. : 1.42922 Max. : 1.534604 Max. : 1.684594
##
## log2_TRFINDEX log2_STRF log2_CTX log2_UI
## Min. :-1.17406 Min. :-1.13215 Min. :-1.48388 Min. :-1.62194
## 1st Qu.:-0.31618 1st Qu.:-0.38289 1st Qu.:-0.26053 1st Qu.:-0.21541
## Median :-0.05137 Median :-0.03349 Median : 0.03195 Median : 0.05194
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.22979 3rd Qu.: 0.27543 3rd Qu.: 0.31093 3rd Qu.: 0.33609
## Max. : 1.74128 Max. : 1.57692 Max. : 1.58523 Max. : 1.06570
##
## log2_CREA log2_HOLOTC log2_HCY log2_MMA
## Min. :-0.8901 Min. :-1.0930 Min. :-1.209351 Min. :-0.9160
## 1st Qu.:-0.4414 1st Qu.:-0.2816 1st Qu.:-0.311324 1st Qu.:-0.3108
## Median :-0.1205 Median : 0.0110 Median :-0.004265 Median :-0.1061
## Mean : 0.0000 Mean : 0.0000 Mean : 0.000000 Mean : 0.0000
## 3rd Qu.: 0.4168 3rd Qu.: 0.2642 3rd Qu.: 0.236266 3rd Qu.: 0.1607
## Max. : 1.3131 Max. : 1.3796 Max. : 1.691373 Max. : 1.9540
##
## log2_VIT_D log2_IGF1 log2_uCa log2_uCCR
## Min. :-1.381056 Min. :-1.0593846 Min. :-0.67077 Min. :-1.14730
## 1st Qu.:-0.366847 1st Qu.:-0.3340024 1st Qu.:-0.45358 1st Qu.:-0.41602
## Median :-0.009946 Median :-0.0003994 Median : 0.03941 Median : 0.03259
## Mean : 0.000000 Mean : 0.0000000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.342680 3rd Qu.: 0.3380598 3rd Qu.: 0.41488 3rd Qu.: 0.37986
## Max. : 1.233721 Max. : 1.0670245 Max. : 1.05893 Max. : 1.09810
##
## log2_uICR
## Min. :-0.99271
## 1st Qu.:-0.33192
## Median :-0.01509
## Mean : 0.00000
## 3rd Qu.: 0.27115
## Max. : 1.59665
## 2.2 Adults
2.2.1 Wrangle data
Open code
for (i in 1:ncol(dat_adult)) {
names(dat_adult)[i] <- str_replace_all(
names(dat_adult)[i], c("^a" = "")
)
}
dat_adult <- dat_adult %>%
rename(
FOLATE = FOLAT,
TIBC = VKFE,
HOLOTC = VIT_AKTB12,
CTX = Cros,
uCCR = Ca_per_Krea,
uICR = I_per_Krea,
uPCR = P_per_Krea,
uCa = Ur_Ca
) %>%
mutate(
log2_FAT = log2(FAT),
log2_FFM = log2(FFM),
log2_TC = log2(TC),
log2_TG = log2(TG),
log2_Se = log2(Se),
log2_Zn = log2(Zn),
log2_TIBC = log2(TIBC),
log2_FERR = log2(FERR),
log2_TRF = log2(TRF),
log2_TRFINDEX = log2(TRFINDEX),
log2_STRF = log2(STRF),
log2_CTX = log2(CTX),
log2_P1NP = log2(P1NP),
log2_UI = log2(UI),
log2_UREA = log2(UREA),
log2_CREA = log2(CREA),
log2_HOLOTC = log2(HOLOTC),
log2_HCY = log2(HCY),
log2_MMA = log2(MMA),
log2_VIT_D = log2(VIT_D),
log2_uCa = log2(uCa),
log2_uCCR = log2(uCCR),
log2_uICR = log2(uICR)
) %>%
mutate(
male = if_else(SEX == "M", 1, 0),
GRP = factor(GRP),
FAM = factor(FAM)
) %>%
select(
## remove irrelevant variables not included to any inference and prediction
-c(SEX, ID),
## remove information about supplementation
-c(SUP_VEG1:SUP_OTH),
## remove non-log2-transformed version when log2(variable) is included
-c(
FAT, FFM, TC, TG, Se, Zn, TIBC, FERR, TRF, TRFINDEX, STRF,
CTX, P1NP, UI, UREA, CREA, HOLOTC, HCY, MMA, VIT_D,
uCa, uCCR, uICR
)
) %>%
select(
FAM, GRP, AGE, male,
everything()
)
summary(dat_adult)
## FAM GRP AGE male BMI
## R1 : 2 OM:50 Min. :21.54 Min. :0.000 Min. :15.89
## R10 : 2 VG:45 1st Qu.:32.47 1st Qu.:0.000 1st Qu.:21.22
## R11 : 2 VN:92 Median :35.32 Median :0.000 Median :23.49
## R12 : 2 Mean :35.54 Mean :0.492 Mean :23.85
## R13 : 2 3rd Qu.:38.27 3rd Qu.:1.000 3rd Qu.:25.98
## R14 : 2 Max. :56.56 Max. :1.000 Max. :34.76
## (Other):175
## WHR HEIGHT MASS HG
## Min. :0.6480 Min. :1.531 Min. : 46.80 Min. :20.75
## 1st Qu.:0.7368 1st Qu.:1.690 1st Qu.: 62.60 1st Qu.:30.75
## Median :0.7843 Median :1.742 Median : 70.40 Median :39.90
## Mean :0.7890 Mean :1.747 Mean : 73.09 Mean :42.12
## 3rd Qu.:0.8341 3rd Qu.:1.804 3rd Qu.: 82.65 3rd Qu.:52.20
## Max. :1.0000 Max. :1.983 Max. :123.50 Max. :76.25
## NA's :4 NA's :1
## PFAT SBP DBP GLY
## Min. : 2.70 Min. : 82.0 Min. : 46.00 Min. :3.100
## 1st Qu.:16.05 1st Qu.:112.5 1st Qu.: 71.00 1st Qu.:4.195
## Median :21.90 Median :123.0 Median : 77.00 Median :4.460
## Mean :22.47 Mean :124.2 Mean : 76.86 Mean :4.454
## 3rd Qu.:28.65 3rd Qu.:135.0 3rd Qu.: 83.00 3rd Qu.:4.720
## Max. :43.90 Max. :179.0 Max. :104.00 Max. :8.410
##
## HDL LDL Ca P
## Min. :0.800 Min. :1.070 Min. :2.180 Min. :0.660
## 1st Qu.:1.215 1st Qu.:2.020 1st Qu.:2.410 1st Qu.:1.005
## Median :1.430 Median :2.440 Median :2.460 Median :1.100
## Mean :1.488 Mean :2.538 Mean :2.456 Mean :1.108
## 3rd Qu.:1.710 3rd Qu.:2.965 3rd Qu.:2.510 3rd Qu.:1.220
## Max. :3.030 Max. :5.960 Max. :2.630 Max. :1.480
##
## Mg FE SATTRF HGB
## Min. :0.6900 Min. : 2.00 Min. : 1.70 Min. : 40.0
## 1st Qu.:0.7700 1st Qu.:14.75 1st Qu.:20.50 1st Qu.:135.5
## Median :0.8100 Median :18.70 Median :27.30 Median :145.0
## Mean :0.8121 Mean :20.04 Mean :29.27 Mean :144.0
## 3rd Qu.:0.8500 3rd Qu.:24.75 3rd Qu.:36.45 3rd Qu.:155.0
## Max. :1.0700 Max. :46.70 Max. :75.20 Max. :182.0
##
## MCV PTH UA FOLATE
## Min. :57.60 Min. :0.800 Min. :138.0 Min. : 3.05
## 1st Qu.:86.00 1st Qu.:2.400 1st Qu.:247.0 1st Qu.: 8.97
## Median :88.40 Median :3.200 Median :296.0 Median :11.59
## Mean :88.61 Mean :3.295 Mean :301.7 Mean :12.46
## 3rd Qu.:91.05 3rd Qu.:4.000 3rd Qu.:349.5 3rd Qu.:15.27
## Max. :99.50 Max. :7.800 Max. :516.0 Max. :24.01
##
## uPCR AL_adult BREAKS log2_FAT
## Min. :0.1565 Min. :0.0000 Min. :0.0000 Min. :0.6781
## 1st Qu.:0.7921 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:3.4983
## Median :1.3387 Median :0.0000 Median :0.0000 Median :3.9260
## Mean :1.4656 Mean :0.2193 Mean :0.4652 Mean :3.9009
## 3rd Qu.:1.9878 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:4.3889
## Max. :4.7073 Max. :1.0000 Max. :1.0000 Max. :5.3255
## NA's :1
## log2_FFM log2_TC log2_TG log2_Se
## Min. :5.289 Min. :1.395 Min. :-1.7859 Min. :-1.47393
## 1st Qu.:5.545 1st Qu.:1.935 1st Qu.:-0.7859 1st Qu.:-0.35845
## Median :5.755 Median :2.138 Median :-0.3585 Median :-0.02915
## Mean :5.791 Mean :2.128 Mean :-0.2534 Mean :-0.08565
## 3rd Qu.:6.060 3rd Qu.:2.328 3rd Qu.: 0.1110 3rd Qu.: 0.23265
## Max. :6.461 Max. :3.064 Max. : 1.8520 Max. : 1.23266
##
## log2_Zn log2_TIBC log2_FERR log2_TRF
## Min. :2.926 Min. :5.728 Min. :-0.5146 Min. :1.084
## 1st Qu.:3.524 1st Qu.:6.000 1st Qu.: 3.9681 1st Qu.:1.350
## Median :3.689 Median :6.109 Median : 4.7970 Median :1.444
## Mean :3.698 Mean :6.125 Mean : 4.8328 Mean :1.470
## 3rd Qu.:3.848 3rd Qu.:6.248 3rd Qu.: 5.5759 3rd Qu.:1.587
## Max. :4.733 Max. :6.907 Max. : 8.6949 Max. :2.257
##
## log2_TRFINDEX log2_STRF log2_CTX log2_P1NP
## Min. :-1.43440 Min. :-0.43440 Min. :-4.3219 Min. :4.433
## 1st Qu.:-0.62149 1st Qu.: 0.06003 1st Qu.:-1.8651 1st Qu.:5.221
## Median :-0.25154 Median : 0.28094 Median :-1.3328 Median :5.558
## Mean :-0.24387 Mean : 0.27996 Mean :-1.3881 Mean :5.643
## 3rd Qu.: 0.08406 3rd Qu.: 0.43028 3rd Qu.:-0.8612 3rd Qu.:6.059
## Max. : 1.62293 Max. : 4.15381 Max. : 0.5945 Max. :7.322
## NA's :2 NA's :1
## log2_UI log2_UREA log2_CREA log2_HOLOTC
## Min. :-0.737 Min. :0.848 Min. :5.322 Min. :4.392
## 1st Qu.: 5.600 1st Qu.:1.766 1st Qu.:5.833 1st Qu.:6.033
## Median : 6.642 Median :2.070 Median :6.022 Median :6.416
## Mean : 6.401 Mean :2.063 Mean :6.036 Mean :6.439
## 3rd Qu.: 7.375 3rd Qu.:2.322 3rd Qu.:6.219 3rd Qu.:6.815
## Max. : 9.647 Max. :4.000 Max. :7.285 Max. :9.313
## NA's :18
## log2_HCY log2_MMA log2_VIT_D log2_uCa
## Min. :2.963 Min. : 6.349 Min. :4.907 Min. :-2.0006
## 1st Qu.:3.620 1st Qu.: 7.178 1st Qu.:5.898 1st Qu.:-1.7251
## Median :3.907 Median : 7.548 Median :6.186 Median :-0.4546
## Mean :3.928 Mean : 7.624 Mean :6.166 Mean :-0.3721
## 3rd Qu.:4.154 3rd Qu.: 7.889 3rd Qu.:6.438 3rd Qu.: 0.4931
## Max. :6.022 Max. :12.613 Max. :7.314 Max. : 3.0548
## NA's :24 NA's :1
## log2_uCCR log2_uICR
## Min. :-6.3997 Min. : 1.481
## 1st Qu.:-3.3701 1st Qu.: 6.223
## Median :-2.7674 Median : 7.170
## Mean :-2.7535 Mean : 7.180
## 3rd Qu.:-2.1385 3rd Qu.: 8.223
## Max. : 0.2224 Max. :11.964
## NA's :1 NA's :192.2.2 Missing vlaues and stadardization
Open code
skimr::skim(dat_adult)| Name | dat_adult |
| Number of rows | 187 |
| Number of columns | 51 |
| _______________________ | |
| Column type frequency: | |
| factor | 2 |
| numeric | 49 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| FAM | 0 | 1 | FALSE | 95 | R1: 2, R10: 2, R11: 2, R12: 2 |
| GRP | 0 | 1 | FALSE | 3 | VN: 92, OM: 50, VG: 45 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| AGE | 0 | 1.00 | 35.54 | 5.35 | 21.54 | 32.47 | 35.32 | 38.27 | 56.56 | ▂▇▇▁▁ |
| male | 0 | 1.00 | 0.49 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| BMI | 0 | 1.00 | 23.85 | 3.59 | 15.89 | 21.22 | 23.49 | 25.98 | 34.76 | ▂▇▇▂▁ |
| WHR | 4 | 0.98 | 0.79 | 0.07 | 0.65 | 0.74 | 0.78 | 0.83 | 1.00 | ▃▇▆▂▁ |
| HEIGHT | 0 | 1.00 | 1.75 | 0.09 | 1.53 | 1.69 | 1.74 | 1.80 | 1.98 | ▂▅▇▃▂ |
| MASS | 0 | 1.00 | 73.09 | 13.97 | 46.80 | 62.60 | 70.40 | 82.65 | 123.50 | ▅▇▆▁▁ |
| HG | 1 | 0.99 | 42.12 | 12.98 | 20.75 | 30.75 | 39.90 | 52.20 | 76.25 | ▇▆▇▅▁ |
| PFAT | 0 | 1.00 | 22.47 | 8.28 | 2.70 | 16.05 | 21.90 | 28.65 | 43.90 | ▁▆▇▅▂ |
| SBP | 0 | 1.00 | 124.20 | 16.90 | 82.00 | 112.50 | 123.00 | 135.00 | 179.00 | ▂▇▇▃▁ |
| DBP | 0 | 1.00 | 76.86 | 9.30 | 46.00 | 71.00 | 77.00 | 83.00 | 104.00 | ▁▂▇▅▁ |
| GLY | 0 | 1.00 | 4.45 | 0.49 | 3.10 | 4.20 | 4.46 | 4.72 | 8.41 | ▂▇▁▁▁ |
| HDL | 0 | 1.00 | 1.49 | 0.40 | 0.80 | 1.21 | 1.43 | 1.71 | 3.03 | ▅▇▃▁▁ |
| LDL | 0 | 1.00 | 2.54 | 0.71 | 1.07 | 2.02 | 2.44 | 2.96 | 5.96 | ▅▇▃▁▁ |
| Ca | 0 | 1.00 | 2.46 | 0.08 | 2.18 | 2.41 | 2.46 | 2.51 | 2.63 | ▁▂▇▇▃ |
| P | 0 | 1.00 | 1.11 | 0.16 | 0.66 | 1.00 | 1.10 | 1.22 | 1.48 | ▁▃▇▆▂ |
| Mg | 0 | 1.00 | 0.81 | 0.07 | 0.69 | 0.77 | 0.81 | 0.85 | 1.07 | ▃▇▃▁▁ |
| FE | 0 | 1.00 | 20.04 | 8.00 | 2.00 | 14.75 | 18.70 | 24.75 | 46.70 | ▂▇▆▂▁ |
| SATTRF | 0 | 1.00 | 29.27 | 12.61 | 1.70 | 20.50 | 27.30 | 36.45 | 75.20 | ▂▇▅▂▁ |
| HGB | 0 | 1.00 | 143.96 | 14.98 | 40.00 | 135.50 | 145.00 | 155.00 | 182.00 | ▁▁▁▇▃ |
| MCV | 0 | 1.00 | 88.61 | 4.24 | 57.60 | 86.00 | 88.40 | 91.05 | 99.50 | ▁▁▁▇▃ |
| PTH | 0 | 1.00 | 3.29 | 1.18 | 0.80 | 2.40 | 3.20 | 4.00 | 7.80 | ▃▇▅▁▁ |
| UA | 0 | 1.00 | 301.66 | 70.78 | 138.00 | 247.00 | 296.00 | 349.50 | 516.00 | ▂▇▇▃▁ |
| FOLATE | 0 | 1.00 | 12.46 | 4.87 | 3.05 | 8.97 | 11.59 | 15.27 | 24.01 | ▂▇▆▂▂ |
| uPCR | 1 | 0.99 | 1.47 | 0.83 | 0.16 | 0.79 | 1.34 | 1.99 | 4.71 | ▇▇▅▁▁ |
| AL_adult | 0 | 1.00 | 0.22 | 0.41 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| BREAKS | 0 | 1.00 | 0.47 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| log2_FAT | 0 | 1.00 | 3.90 | 0.69 | 0.68 | 3.50 | 3.93 | 4.39 | 5.33 | ▁▁▃▇▃ |
| log2_FFM | 0 | 1.00 | 5.79 | 0.29 | 5.29 | 5.55 | 5.75 | 6.06 | 6.46 | ▇▇▅▇▂ |
| log2_TC | 0 | 1.00 | 2.13 | 0.29 | 1.40 | 1.94 | 2.14 | 2.33 | 3.06 | ▂▇▇▃▁ |
| log2_TG | 0 | 1.00 | -0.25 | 0.71 | -1.79 | -0.79 | -0.36 | 0.11 | 1.85 | ▂▇▆▂▁ |
| log2_Se | 0 | 1.00 | -0.09 | 0.49 | -1.47 | -0.36 | -0.03 | 0.23 | 1.23 | ▁▃▇▅▁ |
| log2_Zn | 0 | 1.00 | 3.70 | 0.27 | 2.93 | 3.52 | 3.69 | 3.85 | 4.73 | ▁▆▇▂▁ |
| log2_TIBC | 0 | 1.00 | 6.13 | 0.20 | 5.73 | 6.00 | 6.11 | 6.25 | 6.91 | ▃▇▃▁▁ |
| log2_FERR | 0 | 1.00 | 4.83 | 1.24 | -0.51 | 3.97 | 4.80 | 5.58 | 8.69 | ▁▁▇▆▁ |
| log2_TRF | 0 | 1.00 | 1.47 | 0.20 | 1.08 | 1.35 | 1.44 | 1.59 | 2.26 | ▃▇▃▁▁ |
| log2_TRFINDEX | 2 | 0.99 | -0.24 | 0.53 | -1.43 | -0.62 | -0.25 | 0.08 | 1.62 | ▂▇▆▂▁ |
| log2_STRF | 1 | 0.99 | 0.28 | 0.41 | -0.43 | 0.06 | 0.28 | 0.43 | 4.15 | ▇▂▁▁▁ |
| log2_CTX | 0 | 1.00 | -1.39 | 0.79 | -4.32 | -1.87 | -1.33 | -0.86 | 0.59 | ▁▂▇▇▂ |
| log2_P1NP | 0 | 1.00 | 5.64 | 0.59 | 4.43 | 5.22 | 5.56 | 6.06 | 7.32 | ▃▇▆▃▁ |
| log2_UI | 18 | 0.90 | 6.40 | 1.59 | -0.74 | 5.60 | 6.64 | 7.38 | 9.65 | ▁▁▂▇▂ |
| log2_UREA | 0 | 1.00 | 2.06 | 0.42 | 0.85 | 1.77 | 2.07 | 2.32 | 4.00 | ▁▇▆▁▁ |
| log2_CREA | 0 | 1.00 | 6.04 | 0.29 | 5.32 | 5.83 | 6.02 | 6.22 | 7.29 | ▂▇▆▁▁ |
| log2_HOLOTC | 0 | 1.00 | 6.44 | 0.69 | 4.39 | 6.03 | 6.42 | 6.82 | 9.31 | ▁▇▇▁▁ |
| log2_HCY | 24 | 0.87 | 3.93 | 0.51 | 2.96 | 3.62 | 3.91 | 4.15 | 6.02 | ▃▇▃▁▁ |
| log2_MMA | 0 | 1.00 | 7.62 | 0.68 | 6.35 | 7.18 | 7.55 | 7.89 | 12.61 | ▇▆▁▁▁ |
| log2_VIT_D | 0 | 1.00 | 6.17 | 0.44 | 4.91 | 5.90 | 6.19 | 6.44 | 7.31 | ▁▃▇▆▁ |
| log2_uCa | 1 | 0.99 | -0.37 | 1.31 | -2.00 | -1.73 | -0.45 | 0.49 | 3.05 | ▇▆▅▂▁ |
| log2_uCCR | 1 | 0.99 | -2.75 | 1.01 | -6.40 | -3.37 | -2.77 | -2.14 | 0.22 | ▁▂▇▆▁ |
| log2_uICR | 19 | 0.90 | 7.18 | 1.64 | 1.48 | 6.22 | 7.17 | 8.22 | 11.96 | ▁▂▇▅▁ |
All features has over 3/4 of non-missing values
2.2.2.1 Multiple imputation
Prepare
Open code
dat_adult_toimp <- dat_adult %>% select(-FAM)
init <- mice(dat_adult_toimp, maxit = 0)
## Warning: Number of logged events: 1
init$method
## GRP AGE male BMI WHR
## "" "" "" "" "pmm"
## HEIGHT MASS HG PFAT SBP
## "" "" "pmm" "" ""
## DBP GLY HDL LDL Ca
## "" "" "" "" ""
## P Mg FE SATTRF HGB
## "" "" "" "" ""
## MCV PTH UA FOLATE uPCR
## "" "" "" "" "pmm"
## AL_adult BREAKS log2_FAT log2_FFM log2_TC
## "" "" "" "" ""
## log2_TG log2_Se log2_Zn log2_TIBC log2_FERR
## "" "" "" "" ""
## log2_TRF log2_TRFINDEX log2_STRF log2_CTX log2_P1NP
## "" "pmm" "pmm" "" ""
## log2_UI log2_UREA log2_CREA log2_HOLOTC log2_HCY
## "pmm" "" "" "" "pmm"
## log2_MMA log2_VIT_D log2_uCa log2_uCCR log2_uICR
## "" "" "pmm" "pmm" "pmm"
data_imputed <- run(
expr = mice(
dat_adult_toimp,
method = init$method,
m = 1
),
path = "gitignore/run/data_adult_imputed"
)
dat_adult_imputed <- data.frame(
FAM = dat_adult$FAM,
complete(data_imputed, 'all')[1]$`1`
)2.2.2.2 Rescaling
Open code
dat_adult_imputed <- dat_adult_imputed %>%
mutate(
across(
where(~ is.numeric(.) && length(unique(.)) > 2),
arm::rescale
)
)
summary(dat_adult_imputed)
## FAM GRP AGE male BMI
## R1 : 2 OM:50 Min. :-1.30803 Min. :0.000 Min. :-1.10918
## R10 : 2 VG:45 1st Qu.:-0.28651 1st Qu.:0.000 1st Qu.:-0.36567
## R11 : 2 VN:92 Median :-0.02015 Median :0.000 Median :-0.04986
## R12 : 2 Mean : 0.00000 Mean :0.492 Mean : 0.00000
## R13 : 2 3rd Qu.: 0.25504 3rd Qu.:1.000 3rd Qu.: 0.29673
## R14 : 2 Max. : 1.96479 Max. :1.000 Max. : 1.51954
## (Other):175
## WHR HEIGHT MASS HG
## Min. :-1.02613 Min. :-1.19707 Min. :-0.9408 Min. :-0.8219
## 1st Qu.:-0.37285 1st Qu.:-0.31572 1st Qu.:-0.3753 1st Qu.:-0.4363
## Median :-0.03565 Median :-0.02748 Median :-0.0961 Median :-0.1202
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.32857 3rd Qu.: 0.31619 3rd Qu.: 0.3423 3rd Qu.: 0.3887
## Max. : 1.54683 Max. : 1.30841 Max. : 1.8044 Max. : 1.3178
##
## PFAT SBP DBP GLY
## Min. :-1.19376 Min. :-1.2486 Min. :-1.659628 Min. :-1.368237
## 1st Qu.:-0.38753 1st Qu.:-0.3462 1st Qu.:-0.314955 1st Qu.:-0.262125
## Median :-0.03423 Median :-0.0356 Median : 0.007766 Median : 0.005564
## Mean : 0.00000 Mean : 0.0000 Mean : 0.000000 Mean : 0.000000
## 3rd Qu.: 0.37341 3rd Qu.: 0.3194 3rd Qu.: 0.330487 3rd Qu.: 0.268202
## Max. : 1.29439 Max. : 1.6211 Max. : 1.460013 Max. : 3.995648
##
## HDL LDL Ca P
## Min. :-0.85527 Min. :-1.02864 Min. :-1.7722 Min. :-1.415
## 1st Qu.:-0.33963 1st Qu.:-0.36317 1st Qu.:-0.2950 1st Qu.:-0.325
## Median :-0.07249 Median :-0.06896 Median : 0.0261 Median :-0.025
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000 Mean : 0.000
## 3rd Qu.: 0.27541 3rd Qu.: 0.29880 3rd Qu.: 0.3472 3rd Qu.: 0.354
## Max. : 1.91552 Max. : 2.39677 Max. : 1.1179 Max. : 1.175
##
## Mg FE SATTRF HGB
## Min. :-0.91127 Min. :-1.12804 Min. :-1.09328 Min. :-3.4692
## 1st Qu.:-0.31440 1st Qu.:-0.33095 1st Qu.:-0.34777 1st Qu.:-0.2822
## Median :-0.01596 Median :-0.08401 Median :-0.07812 Median : 0.0348
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.28248 3rd Qu.: 0.29421 3rd Qu.: 0.28472 3rd Qu.: 0.3685
## Max. : 1.92389 Max. : 1.66645 Max. : 1.82133 Max. : 1.2696
##
## MCV PTH UA FOLATE
## Min. :-3.65974 Min. :-1.06081 Min. :-1.15606 Min. :-0.96640
## 1st Qu.:-0.30764 1st Qu.:-0.38044 1st Qu.:-0.38610 1st Qu.:-0.35870
## Median :-0.02436 Median :-0.04025 Median :-0.03997 Median :-0.08975
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.28842 3rd Qu.: 0.29994 3rd Qu.: 0.33795 3rd Qu.: 0.28749
## Max. : 1.28579 Max. : 1.91583 Max. : 1.51409 Max. : 1.18518
##
## uPCR AL_adult BREAKS log2_FAT
## Min. :-0.79306 Min. :0.0000 Min. :0.0000 Min. :-2.3394
## 1st Qu.:-0.40597 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:-0.2923
## Median :-0.08154 Median :0.0000 Median :0.0000 Median : 0.0182
## Mean : 0.00000 Mean :0.2193 Mean :0.4652 Mean : 0.0000
## 3rd Qu.: 0.31714 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.: 0.3542
## Max. : 1.96700 Max. :1.0000 Max. :1.0000 Max. : 1.0341
##
## log2_FFM log2_TC log2_TG log2_Se
## Min. :-0.85542 Min. :-1.27470 Min. :-1.07569 Min. :-1.42194
## 1st Qu.:-0.41825 1st Qu.:-0.33509 1st Qu.:-0.37376 1st Qu.:-0.27942
## Median :-0.06078 Median : 0.01621 Median :-0.07373 Median : 0.05788
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.45908 3rd Qu.: 0.34690 3rd Qu.: 0.25581 3rd Qu.: 0.32602
## Max. : 1.14396 Max. : 1.62630 Max. : 1.47786 Max. : 1.35028
##
## log2_Zn log2_TIBC log2_FERR log2_TRF
## Min. :-1.41578 Min. :-0.99967 Min. :-2.15715 Min. :-0.97211
## 1st Qu.:-0.31993 1st Qu.:-0.31549 1st Qu.:-0.34884 1st Qu.:-0.30136
## Median :-0.01599 Median :-0.04259 Median :-0.01445 Median :-0.06695
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.27504 3rd Qu.: 0.30795 3rd Qu.: 0.29976 3rd Qu.: 0.29496
## Max. : 1.89867 Max. : 1.96499 Max. : 1.55796 Max. : 1.98082
##
## log2_TRFINDEX log2_STRF log2_CTX log2_P1NP
## Min. :-1.10778 Min. :-0.877703 Min. :-1.8463 Min. :-1.02184
## 1st Qu.:-0.34728 1st Qu.:-0.266816 1st Qu.:-0.3002 1st Qu.:-0.35636
## Median :-0.01598 Median : 0.007279 Median : 0.0348 Median :-0.07218
## Mean : 0.00000 Mean : 0.000000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.29378 3rd Qu.: 0.186510 3rd Qu.: 0.3316 3rd Qu.: 0.35073
## Max. : 1.71417 Max. : 4.751818 Max. : 1.2477 Max. : 1.41750
##
## log2_UI log2_UREA log2_CREA log2_HOLOTC
## Min. :-2.24120 Min. :-1.432231 Min. :-1.21223 Min. :-1.47624
## 1st Qu.:-0.25347 1st Qu.:-0.350970 1st Qu.:-0.34506 1st Qu.:-0.29262
## Median : 0.09545 Median : 0.008282 Median :-0.02349 Median :-0.01657
## Mean : 0.00000 Mean : 0.000000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.29693 3rd Qu.: 0.304705 3rd Qu.: 0.31045 3rd Qu.: 0.27113
## Max. : 0.99762 Max. : 2.282208 Max. : 2.12003 Max. : 2.07264
##
## log2_HCY log2_MMA log2_VIT_D log2_uCa
## Min. :-0.958835 Min. :-0.94318 Min. :-1.43672 Min. :-0.62121
## 1st Qu.:-0.307167 1st Qu.:-0.32979 1st Qu.:-0.30517 1st Qu.:-0.51130
## Median :-0.004007 Median :-0.05622 Median : 0.02283 Median :-0.04482
## Mean : 0.000000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.249839 3rd Qu.: 0.19536 3rd Qu.: 0.31051 3rd Qu.: 0.33192
## Max. : 2.075618 Max. : 3.68813 Max. : 1.30999 Max. : 1.31305
##
## log2_uCCR log2_uICR
## Min. :-1.804807 Min. :-1.78656
## 1st Qu.:-0.304338 1st Qu.:-0.30019
## Median :-0.006875 Median :-0.01302
## Mean : 0.000000 Mean : 0.00000
## 3rd Qu.: 0.303717 3rd Qu.: 0.33061
## Max. : 1.472229 Max. : 1.47788
## 3 Diet prediction
(OM, VN, VG determination)
All prediction models will account for the variable FAM as the grouping (block) variable, ensuring that all subject of a single family will be part of the same fold during each bootstrap iteration.
3.1 All children
3.1.1 Baseline model
Following models will be trained to discriminate between two diets, based solely on log2_age, male, BreastFeed_total_stopped, BreastFeed_full_stopped and interaction BreastFeed_full_duration.
3.1.1.1 Prepare data
Open code
dat_pred <- dat_child_all_imputed %>%
select(GRP, FAM, log2_age, male, BreastFeed_full_duration,
BreastFeed_full_stopped, BreastFeed_total_stopped, BiW)
nfam <- dat_pred$FAM %>% factor() %>% levels() %>% length()
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "log2_age" "male"
## [3] "BreastFeed_full_duration" "BreastFeed_full_stopped"
## [5] "BreastFeed_total_stopped" "BiW"
dat_outcome <- dat_pred %>% select(GRP) %>% as.matrix3.1.1.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.5809524
## Fit and validate model
modelac <- 'pred_baseline_child_all_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.4039723
## auc 0.6350596
## auc_optimism_corrected 0.5372320
## auc_optimism_corrected_CIL 0.3586458
## auc_optimism_corrected_CIU 0.6872924
## accuracy 0.5809524
## accuracy_optimism_corrected 0.5327153
## accuracy_optimism_corrected_CIL 0.3468919
## accuracy_optimism_corrected_CIU 0.6775735
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism
## Min. :0.5000 Min. :0.2742 Min. :-0.134138
## 1st Qu.:0.6168 1st Qu.:0.5000 1st Qu.: 0.001095
## Median :0.6739 Median :0.5284 Median : 0.108890
## Mean :0.6582 Mean :0.5372 Mean : 0.120975
## 3rd Qu.:0.7238 3rd Qu.:0.5890 3rd Qu.: 0.202375
## Max. :0.8917 Max. :0.7717 Max. : 0.527510
## accuracy_resamp_test accuracy_resamp_validation accuracy_optimism
## Min. :0.5049 Min. :0.2000 Min. :-0.22660
## 1st Qu.:0.5941 1st Qu.:0.4773 1st Qu.: 0.01389
## Median :0.6281 Median :0.5366 Median : 0.09871
## Mean :0.6304 Mean :0.5327 Mean : 0.09765
## 3rd Qu.:0.6667 3rd Qu.:0.5952 3rd Qu.: 0.17749
## Max. :0.8317 Max. :0.7576 Max. : 0.44119
## Model coefficients
get(modelac)$betas
## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.06622559
## male .
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped -0.07857001
## BiW .3.1.1.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4567901
## Fit and validate model
modelac <- 'pred_baseline_child_all_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.0000000
## lambda 0.1504131
## auc 0.7542998
## auc_optimism_corrected 0.5749206
## auc_optimism_corrected_CIL 0.3285965
## auc_optimism_corrected_CIU 0.7837943
## accuracy 0.6790123
## accuracy_optimism_corrected 0.5097507
## accuracy_optimism_corrected_CIL 0.2815476
## accuracy_optimism_corrected_CIU 0.7030882
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.1732 Min. :-0.1778 Min. :0.5065
## 1st Qu.:0.7422 1st Qu.:0.5000 1st Qu.: 0.1088 1st Qu.:0.6657
## Median :0.7898 Median :0.5825 Median : 0.1972 Median :0.7089
## Mean :0.7718 Mean :0.5749 Mean : 0.1968 Mean :0.7024
## 3rd Qu.:0.8227 3rd Qu.:0.6597 3rd Qu.: 0.2802 3rd Qu.:0.7468
## Max. :0.9410 Max. :0.8681 Max. : 0.5914 Max. :0.8831
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1765 Min. :-0.1658
## 1st Qu.:0.4409 1st Qu.: 0.1102
## Median :0.5172 Median : 0.1903
## Mean :0.5098 Mean : 0.1927
## 3rd Qu.:0.5862 3rd Qu.: 0.2778
## Max. :0.8182 Max. : 0.5610
## Model coefficients
get(modelac)$betas
## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age 0.39328720
## male 0.19015857
## BreastFeed_full_duration 0.71146159
## BreastFeed_full_stopped -1.05623482
## BreastFeed_total_stopped -0.58807096
## BiW -0.067137963.1.1.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.622449
## Fit and validate model
modelac <- 'pred_baseline_child_all_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.8000000
## lambda 0.1067586
## auc 0.6637129
## auc_optimism_corrected 0.5923943
## auc_optimism_corrected_CIL 0.3866566
## auc_optimism_corrected_CIU 0.7576386
## accuracy 0.6224490
## accuracy_optimism_corrected 0.6002464
## accuracy_optimism_corrected_CIL 0.4475267
## accuracy_optimism_corrected_CIU 0.7567568
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.1781 Min. :-0.15986 Min. :0.2400
## 1st Qu.:0.6841 1st Qu.:0.5235 1st Qu.: 0.03207 1st Qu.:0.6504
## Median :0.7391 Median :0.6029 Median : 0.12478 Median :0.6848
## Mean :0.7221 Mean :0.5924 Mean : 0.12972 Mean :0.6866
## 3rd Qu.:0.7848 3rd Qu.:0.6619 3rd Qu.: 0.20828 3rd Qu.:0.7263
## Max. :0.8872 Max. :0.8733 Max. : 0.55336 Max. :0.8416
## accuracy_resamp_validation accuracy_optimism
## Min. :0.3889 Min. :-0.29846
## 1st Qu.:0.5526 1st Qu.: 0.01585
## Median :0.6000 Median : 0.08727
## Mean :0.6002 Mean : 0.08634
## 3rd Qu.:0.6471 3rd Qu.: 0.16024
## Max. :0.8667 Max. : 0.41571
## Model coefficients
get(modelac)$betas
## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.29889
## male .
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped .
## BiW .3.1.2 Complete model
Baseline + clinical characteristics
3.1.2.1 Prepare data
Open code
dat_pred <- dat_child_all_imputed %>%
select(GRP, FAM, log2_age, male, BreastFeed_full_duration,
BreastFeed_full_stopped, BreastFeed_total_stopped,
BiW,
MASS_Perc:log2_uICR)
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "log2_age" "male"
## [3] "BreastFeed_full_duration" "BreastFeed_full_stopped"
## [5] "BreastFeed_total_stopped" "BiW"
## [7] "MASS_Perc" "HEIGHT_Perc"
## [9] "BMI_PERC" "M_per_H_PERC"
## [11] "GLY" "TC"
## [13] "HDL" "LDL"
## [15] "Ca" "P"
## [17] "Mg" "FE"
## [19] "TIBC" "TRF"
## [21] "SATTRF" "HGB"
## [23] "MCV" "PTH"
## [25] "P1NP" "UREA"
## [27] "UA" "FOLATE"
## [29] "uPCR" "AL_child"
## [31] "log2_TG" "log2_Se"
## [33] "log2_Zn" "log2_FERR"
## [35] "log2_TRFINDEX" "log2_STRF"
## [37] "log2_CTX" "log2_UI"
## [39] "log2_CREA" "log2_HOLOTC"
## [41] "log2_HCY" "log2_MMA"
## [43] "log2_VIT_D" "log2_IGF1"
## [45] "log2_uCa" "log2_uCCR"
## [47] "log2_uICR"3.1.2.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.5809524
## Fit and validate model
modelac <- 'pred_complet_child_all_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
get(modelac)$plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Open code
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.0999577
## auc 0.9828614
## auc_optimism_corrected 0.8646881
## auc_optimism_corrected_CIL 0.7430309
## auc_optimism_corrected_CIU 0.9522917
## accuracy 0.9333333
## accuracy_optimism_corrected 0.7781709
## accuracy_optimism_corrected_CIL 0.6363991
## accuracy_optimism_corrected_CIU 0.9004756
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.9551 Min. :0.6565 Min. :-0.01226 Min. :0.8796
## 1st Qu.:0.9943 1st Qu.:0.8333 1st Qu.: 0.09301 1st Qu.:0.9623
## Median :0.9990 Median :0.8721 Median : 0.12419 Median :0.9811
## Mean :0.9957 Mean :0.8647 Mean : 0.13099 Mean :0.9756
## 3rd Qu.:1.0000 3rd Qu.:0.9033 3rd Qu.: 0.16485 3rd Qu.:1.0000
## Max. :1.0000 Max. :0.9931 Max. : 0.34345 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.5208 Min. :-0.08761
## 1st Qu.:0.7368 1st Qu.: 0.14786
## Median :0.7838 Median : 0.19982
## Mean :0.7782 Mean : 0.19743
## 3rd Qu.:0.8235 3rd Qu.: 0.24390
## Max. :0.9722 Max. : 0.47917
## Model coefficients
get(modelac)$betas
## 47 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.023259913
## male 0.034639697
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped -0.187358974
## BiW -0.225673112
## MASS_Perc -0.124980706
## HEIGHT_Perc -0.352304701
## BMI_PERC .
## M_per_H_PERC .
## GLY 0.008439785
## TC .
## HDL 0.037182485
## LDL -0.146417229
## Ca .
## P .
## Mg 0.271171726
## FE .
## TIBC -0.213925698
## TRF -0.321646059
## SATTRF 0.116861635
## HGB -0.139453747
## MCV 0.652158784
## PTH 0.129721503
## P1NP 0.081592097
## UREA -0.258241015
## UA -0.029756305
## FOLATE 0.856689101
## uPCR -1.103703323
## AL_child -0.140009893
## log2_TG .
## log2_Se -0.027280307
## log2_Zn -0.155872826
## log2_FERR -0.665202312
## log2_TRFINDEX 0.250201535
## log2_STRF .
## log2_CTX .
## log2_UI -0.084845177
## log2_CREA -0.372434731
## log2_HOLOTC 0.484119828
## log2_HCY .
## log2_MMA -0.953253014
## log2_VIT_D 0.392936461
## log2_IGF1 -0.355548896
## log2_uCa -0.023014066
## log2_uCCR .
## log2_uICR .3.1.2.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4567901
## Fit and validate model
modelac <- 'pred_complet_child_all_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.6000000
## lambda 0.1728403
## auc 0.7979115
## auc_optimism_corrected 0.5908125
## auc_optimism_corrected_CIL 0.3239306
## auc_optimism_corrected_CIU 0.8460016
## accuracy 0.7160494
## accuracy_optimism_corrected 0.5438350
## accuracy_optimism_corrected_CIL 0.3214286
## accuracy_optimism_corrected_CIU 0.7500000
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.2173 Min. :-0.1640 Min. :0.5658
## 1st Qu.:0.9085 1st Qu.:0.5048 1st Qu.: 0.2347 1st Qu.:0.8292
## Median :0.9500 Median :0.5887 Median : 0.3690 Median :0.8902
## Mean :0.9385 Mean :0.5908 Mean : 0.3477 Mean :0.8790
## 3rd Qu.:0.9845 3rd Qu.:0.6815 3rd Qu.: 0.4635 3rd Qu.:0.9452
## Max. :1.0000 Max. :0.9667 Max. : 0.7811 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1429 Min. :-0.1357
## 1st Qu.:0.4682 1st Qu.: 0.2278
## Median :0.5484 Median : 0.3466
## Mean :0.5438 Mean : 0.3352
## 3rd Qu.:0.6167 3rd Qu.: 0.4450
## Max. :0.8750 Max. : 0.8275
## Model coefficients
get(modelac)$betas
## 47 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age .
## male .
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped .
## BiW .
## MASS_Perc .
## HEIGHT_Perc .
## BMI_PERC .
## M_per_H_PERC .
## GLY .
## TC .
## HDL .
## LDL .
## Ca .
## P .
## Mg .
## FE .
## TIBC .
## TRF .
## SATTRF .
## HGB .
## MCV .
## PTH .
## P1NP .
## UREA .
## UA .
## FOLATE 0.90209325
## uPCR -0.04462848
## AL_child .
## log2_TG .
## log2_Se .
## log2_Zn .
## log2_FERR .
## log2_TRFINDEX .
## log2_STRF .
## log2_CTX .
## log2_UI .
## log2_CREA .
## log2_HOLOTC .
## log2_HCY .
## log2_MMA .
## log2_VIT_D .
## log2_IGF1 .
## log2_uCa .
## log2_uCCR .
## log2_uICR .3.1.2.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.622449
## Fit and validate model
modelac <- 'pred_complet_child_all_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 1.00000000
## lambda 0.05240856
## auc 0.92822331
## auc_optimism_corrected 0.70426316
## auc_optimism_corrected_CIL 0.52312653
## auc_optimism_corrected_CIU 0.87598815
## accuracy 0.80612245
## accuracy_optimism_corrected 0.66692062
## accuracy_optimism_corrected_CIL 0.50000000
## accuracy_optimism_corrected_CIU 0.81422697
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism
## Min. :0.7897 Min. :0.3720 Min. :-0.002957
## 1st Qu.:0.9763 1st Qu.:0.6477 1st Qu.: 0.221844
## Median :0.9915 Median :0.7097 Median : 0.275220
## Mean :0.9836 Mean :0.7043 Mean : 0.279333
## 3rd Qu.:0.9991 3rd Qu.:0.7622 3rd Qu.: 0.340465
## Max. :1.0000 Max. :0.9592 Max. : 0.607497
## accuracy_resamp_test accuracy_resamp_validation accuracy_optimism
## Min. :0.6383 Min. :0.3889 Min. :-0.04508
## 1st Qu.:0.9164 1st Qu.:0.6130 1st Qu.: 0.21897
## Median :0.9500 Median :0.6667 Median : 0.27420
## Mean :0.9419 Mean :0.6669 Mean : 0.27501
## 3rd Qu.:0.9796 3rd Qu.:0.7227 3rd Qu.: 0.33346
## Max. :1.0000 Max. :0.9286 Max. : 0.54862
## Model coefficients
get(modelac)$betas
## 47 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.37604128
## male .
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped .
## BiW .
## MASS_Perc .
## HEIGHT_Perc -0.42563039
## BMI_PERC .
## M_per_H_PERC .
## GLY .
## TC .
## HDL 0.02568811
## LDL -0.01840456
## Ca .
## P .
## Mg .
## FE .
## TIBC .
## TRF .
## SATTRF .
## HGB .
## MCV 0.52833146
## PTH 0.47085025
## P1NP .
## UREA .
## UA .
## FOLATE .
## uPCR -0.79802410
## AL_child .
## log2_TG 0.02731084
## log2_Se .
## log2_Zn .
## log2_FERR -0.06851785
## log2_TRFINDEX 0.30304989
## log2_STRF .
## log2_CTX .
## log2_UI .
## log2_CREA .
## log2_HOLOTC 0.75759266
## log2_HCY .
## log2_MMA -0.49930782
## log2_VIT_D .
## log2_IGF1 -1.00932674
## log2_uCa .
## log2_uCCR .
## log2_uICR .3.1.3 Reduced
Baseline + clinical characteristics not affected with vitamins supplementation
3.1.3.1 Prepare data
Open code
dat_pred <- dat_child_all_imputed %>%
select(GRP, FAM, log2_age, male, BreastFeed_full_duration,
BreastFeed_full_stopped, BreastFeed_total_stopped,
BiW,
MASS_Perc:Mg,
log2_TG:log2_Zn,
HGB:UA,
log2_CTX, log2_CREA, log2_IGF1
)
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "log2_age" "male"
## [3] "BreastFeed_full_duration" "BreastFeed_full_stopped"
## [5] "BreastFeed_total_stopped" "BiW"
## [7] "MASS_Perc" "HEIGHT_Perc"
## [9] "BMI_PERC" "M_per_H_PERC"
## [11] "GLY" "TC"
## [13] "HDL" "LDL"
## [15] "Ca" "P"
## [17] "Mg" "log2_TG"
## [19] "log2_Se" "log2_Zn"
## [21] "HGB" "MCV"
## [23] "PTH" "P1NP"
## [25] "UREA" "UA"
## [27] "log2_CTX" "log2_CREA"
## [29] "log2_IGF1"3.1.3.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.5809524
## Fit and validate model
modelac <- 'pred_complex_child_all_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
get(modelac)$plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Open code
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.0000000
## lambda 0.6809560
## auc 0.8386736
## auc_optimism_corrected 0.6525305
## auc_optimism_corrected_CIL 0.4571037
## auc_optimism_corrected_CIU 0.7993206
## accuracy 0.7809524
## accuracy_optimism_corrected 0.6243409
## accuracy_optimism_corrected_CIL 0.4534091
## accuracy_optimism_corrected_CIU 0.7692308
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.2737 Min. :-0.07665 Min. :0.5701
## 1st Qu.:0.8850 1st Qu.:0.5949 1st Qu.: 0.17870 1st Qu.:0.8076
## Median :0.9131 Median :0.6616 Median : 0.24716 Median :0.8500
## Mean :0.9089 Mean :0.6525 Mean : 0.25639 Mean :0.8454
## 3rd Qu.:0.9421 3rd Qu.:0.7162 3rd Qu.: 0.33753 3rd Qu.:0.8879
## Max. :1.0000 Max. :0.9261 Max. : 0.64381 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.3500 Min. :-0.1079
## 1st Qu.:0.5714 1st Qu.: 0.1400
## Median :0.6250 Median : 0.2216
## Mean :0.6243 Mean : 0.2210
## 3rd Qu.:0.6832 3rd Qu.: 0.2973
## Max. :0.8649 Max. : 0.5675
## Model coefficients
get(modelac)$betas
## 29 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.110371269
## male 0.053478581
## BreastFeed_full_duration 0.106491593
## BreastFeed_full_stopped -0.056983757
## BreastFeed_total_stopped -0.130809954
## BiW -0.062498656
## MASS_Perc -0.137379843
## HEIGHT_Perc -0.118087003
## BMI_PERC -0.044998977
## M_per_H_PERC -0.049822380
## GLY 0.010268256
## TC -0.111735304
## HDL 0.025674988
## LDL -0.152029070
## Ca -0.008305915
## P -0.082862856
## Mg 0.201323185
## log2_TG 0.022002611
## log2_Se -0.045072287
## log2_Zn -0.132387542
## HGB -0.084029793
## MCV 0.283508908
## PTH 0.081431561
## P1NP 0.138678669
## UREA -0.289097234
## UA -0.110530244
## log2_CTX 0.097580361
## log2_CREA -0.179034229
## log2_IGF1 -0.2222600633.1.3.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4567901
## Fit and validate model
modelac <- 'pred_complex_child_all_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.6000000
## lambda 0.1791455
## auc 0.5000000
## auc_optimism_corrected 0.4245962
## auc_optimism_corrected_CIL 0.2087214
## auc_optimism_corrected_CIU 0.6402794
## accuracy 0.5432099
## accuracy_optimism_corrected 0.4229261
## accuracy_optimism_corrected_CIL 0.2366246
## accuracy_optimism_corrected_CIU 0.5956818
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.1324 Min. :-0.04312 Min. :0.5057
## 1st Qu.:0.8268 1st Qu.:0.3550 1st Qu.: 0.36449 1st Qu.:0.7124
## Median :0.8900 Median :0.4269 Median : 0.46237 Median :0.7975
## Mean :0.8636 Mean :0.4246 Mean : 0.43898 Mean :0.7823
## 3rd Qu.:0.9376 3rd Qu.:0.5000 3rd Qu.: 0.54385 3rd Qu.:0.8593
## Max. :1.0000 Max. :0.7273 Max. : 0.79429 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1613 Min. :-0.1658
## 1st Qu.:0.3585 1st Qu.: 0.2738
## Median :0.4202 Median : 0.3767
## Mean :0.4229 Mean : 0.3594
## 3rd Qu.:0.4839 3rd Qu.: 0.4559
## Max. :0.6786 Max. : 0.7059
## Model coefficients
get(modelac)$betas
## 29 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age 0
## male .
## BreastFeed_full_duration .
## BreastFeed_full_stopped .
## BreastFeed_total_stopped .
## BiW .
## MASS_Perc .
## HEIGHT_Perc .
## BMI_PERC .
## M_per_H_PERC .
## GLY .
## TC .
## HDL .
## LDL .
## Ca .
## P .
## Mg .
## log2_TG .
## log2_Se .
## log2_Zn .
## HGB .
## MCV .
## PTH .
## P1NP .
## UREA .
## UA .
## log2_CTX .
## log2_CREA .
## log2_IGF1 .3.1.3.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.622449
## Fit and validate model
modelac <- 'pred_complex_child_all_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.0000000
## lambda 0.7235655
## auc 0.8223305
## auc_optimism_corrected 0.6437341
## auc_optimism_corrected_CIL 0.4675493
## auc_optimism_corrected_CIU 0.7907643
## accuracy 0.7142857
## accuracy_optimism_corrected 0.6341987
## accuracy_optimism_corrected_CIL 0.4736842
## accuracy_optimism_corrected_CIU 0.7765848
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.3273 Min. :-0.0257 Min. :0.6383
## 1st Qu.:0.8980 1st Qu.:0.5858 1st Qu.: 0.2044 1st Qu.:0.8122
## Median :0.9283 Median :0.6504 Median : 0.2744 Median :0.8601
## Mean :0.9211 Mean :0.6437 Mean : 0.2774 Mean :0.8538
## 3rd Qu.:0.9571 3rd Qu.:0.7053 3rd Qu.: 0.3508 3rd Qu.:0.9000
## Max. :1.0000 Max. :0.8878 Max. : 0.5962 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.3889 Min. :-0.06015
## 1st Qu.:0.5833 1st Qu.: 0.14928
## Median :0.6389 Median : 0.22152
## Mean :0.6342 Mean : 0.21961
## 3rd Qu.:0.6897 3rd Qu.: 0.29115
## Max. :0.8333 Max. : 0.52632
## Model coefficients
get(modelac)$betas
## 29 x 1 sparse Matrix of class "dgCMatrix"
## s0
## log2_age -0.150337798
## male -0.015480652
## BreastFeed_full_duration -0.040635603
## BreastFeed_full_stopped 0.205739460
## BreastFeed_total_stopped 0.009038895
## BiW -0.031318165
## MASS_Perc -0.091286922
## HEIGHT_Perc -0.162930580
## BMI_PERC 0.012574086
## M_per_H_PERC 0.023907659
## GLY 0.012894483
## TC -0.046958074
## HDL 0.116489341
## LDL -0.205738166
## Ca 0.025643757
## P 0.018195206
## Mg 0.011956604
## log2_TG 0.234485157
## log2_Se -0.020382580
## log2_Zn -0.124665155
## HGB -0.116479185
## MCV 0.240224485
## PTH 0.215368558
## P1NP 0.100360490
## UREA -0.141122560
## UA -0.036142120
## log2_CTX 0.084311132
## log2_CREA -0.148321250
## log2_IGF1 -0.2179917793.2 Adults
3.2.1 Baseline model
Following models will be trained to discriminate between two diets, based solely on AGE and SEX.
3.2.1.1 Prepare data
Open code
dat_pred <- dat_adult_imputed %>%
select(GRP, FAM, AGE, male)
nfam <- dat_pred$FAM %>% factor() %>% levels() %>% length()
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "AGE" "male"
dat_outcome <- dat_pred %>% select(GRP) %>% as.matrix3.2.1.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6478873
## Fit and validate model
modelac <- 'pred_baseline_adult_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
get(modelac)$plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Open code
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.4340878
## auc 0.5000000
## auc_optimism_corrected 0.5423400
## auc_optimism_corrected_CIL 0.3829571
## auc_optimism_corrected_CIU 0.7098909
## accuracy 0.6478873
## accuracy_optimism_corrected 0.6371038
## accuracy_optimism_corrected_CIL 0.4797907
## accuracy_optimism_corrected_CIU 0.7724352
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.4772 Min. :0.2129 Min. :-0.23936 Min. :0.5035
## 1st Qu.:0.5000 1st Qu.:0.5000 1st Qu.: 0.00000 1st Qu.:0.6224
## Median :0.5587 Median :0.5021 Median : 0.00000 Median :0.6572
## Mean :0.5774 Mean :0.5423 Mean : 0.03508 Mean :0.6578
## 3rd Qu.:0.6479 3rd Qu.:0.5990 3rd Qu.: 0.09408 3rd Qu.:0.6901
## Max. :0.7760 Max. :0.7733 Max. : 0.38934 Max. :0.8440
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1695 Min. :-0.39192
## 1st Qu.:0.5889 1st Qu.:-0.06042
## Median :0.6400 Median : 0.01420
## Mean :0.6371 Mean : 0.02068
## 3rd Qu.:0.6897 3rd Qu.: 0.10019
## Max. :0.9130 Max. : 0.43414
## Model coefficients
get(modelac)$betas
## 2 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE 0
## male .3.2.1.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4736842
## Fit and validate model
modelac <- 'pred_baseline_adult_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.60000000
## lambda 0.01329009
## auc 0.52333333
## auc_optimism_corrected 0.46703618
## auc_optimism_corrected_CIL 0.30289855
## auc_optimism_corrected_CIU 0.58972902
## accuracy 0.52631579
## accuracy_optimism_corrected 0.45329487
## accuracy_optimism_corrected_CIL 0.26000000
## accuracy_optimism_corrected_CIU 0.62043269
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.4300 Min. :0.1683 Min. :-0.20455 Min. :0.4167
## 1st Qu.:0.5000 1st Qu.:0.4294 1st Qu.: 0.00000 1st Qu.:0.5263
## Median :0.5137 Median :0.5000 Median : 0.01497 Median :0.5579
## Mean :0.5332 Mean :0.4670 Mean : 0.06621 Mean :0.5654
## 3rd Qu.:0.5586 3rd Qu.:0.5000 3rd Qu.: 0.12488 3rd Qu.:0.5895
## Max. :0.7344 Max. :0.7151 Max. : 0.52763 Max. :0.7419
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1515 Min. :-0.19443
## 1st Qu.:0.3889 1st Qu.: 0.02412
## Median :0.4615 Median : 0.09530
## Mean :0.4533 Mean : 0.11213
## 3rd Qu.:0.5238 3rd Qu.: 0.20063
## Max. :0.6897 Max. : 0.55682
## Model coefficients
get(modelac)$betas
## 2 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE -4.12632e-17
## male .3.2.1.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6715328
## Fit and validate model
modelac <- 'pred_baseline_adult_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 1.00000000
## lambda 0.08085931
## auc 0.50000000
## auc_optimism_corrected 0.54883991
## auc_optimism_corrected_CIL 0.40903409
## auc_optimism_corrected_CIU 0.73334865
## accuracy 0.67153285
## accuracy_optimism_corrected 0.66674384
## accuracy_optimism_corrected_CIL 0.52882908
## accuracy_optimism_corrected_CIU 0.80701754
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.4699 Min. :0.1779 Min. :-0.29205 Min. :0.5036
## 1st Qu.:0.5000 1st Qu.:0.5000 1st Qu.: 0.00000 1st Qu.:0.6350
## Median :0.5385 Median :0.5167 Median : 0.00000 Median :0.6752
## Mean :0.5743 Mean :0.5488 Mean : 0.02549 Mean :0.6728
## 3rd Qu.:0.6488 3rd Qu.:0.5996 3rd Qu.: 0.07911 3rd Qu.:0.7122
## Max. :0.7632 Max. :0.8387 Max. : 0.37786 Max. :0.8406
## accuracy_resamp_validation accuracy_optimism
## Min. :0.4048 Min. :-0.353546
## 1st Qu.:0.6205 1st Qu.:-0.070669
## Median :0.6604 Median : 0.012268
## Mean :0.6667 Mean : 0.006068
## 3rd Qu.:0.7143 3rd Qu.: 0.085948
## Max. :0.8679 Max. : 0.435818
## Model coefficients
get(modelac)$betas
## 2 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE 0
## male .3.2.2 Complete model
Baseline + clinical characteristics
3.2.2.1 Prepare data
Open code
dat_pred <- dat_adult_imputed
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "AGE" "male" "BMI" "WHR"
## [5] "HEIGHT" "MASS" "HG" "PFAT"
## [9] "SBP" "DBP" "GLY" "HDL"
## [13] "LDL" "Ca" "P" "Mg"
## [17] "FE" "SATTRF" "HGB" "MCV"
## [21] "PTH" "UA" "FOLATE" "uPCR"
## [25] "AL_adult" "BREAKS" "log2_FAT" "log2_FFM"
## [29] "log2_TC" "log2_TG" "log2_Se" "log2_Zn"
## [33] "log2_TIBC" "log2_FERR" "log2_TRF" "log2_TRFINDEX"
## [37] "log2_STRF" "log2_CTX" "log2_P1NP" "log2_UI"
## [41] "log2_UREA" "log2_CREA" "log2_HOLOTC" "log2_HCY"
## [45] "log2_MMA" "log2_VIT_D" "log2_uCa" "log2_uCCR"
## [49] "log2_uICR"3.2.2.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6478873
## Fit and validate model
modelac <- 'pred_complet_adult_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
get(modelac)$plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Open code
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.1200082
## auc 0.9593478
## auc_optimism_corrected 0.8547846
## auc_optimism_corrected_CIL 0.7490948
## auc_optimism_corrected_CIU 0.9429849
## accuracy 0.8943662
## accuracy_optimism_corrected 0.7744900
## accuracy_optimism_corrected_CIL 0.6525954
## accuracy_optimism_corrected_CIU 0.8800000
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism
## Min. :0.9347 Min. :0.6540 Min. :-0.006052
## 1st Qu.:0.9844 1st Qu.:0.8242 1st Qu.: 0.097296
## Median :0.9904 Median :0.8561 Median : 0.132035
## Mean :0.9890 Mean :0.8548 Mean : 0.134212
## 3rd Qu.:0.9960 3rd Qu.:0.8906 3rd Qu.: 0.170560
## Max. :1.0000 Max. :0.9768 Max. : 0.336612
## accuracy_resamp_test accuracy_resamp_validation accuracy_optimism
## Min. :0.8462 Min. :0.5690 Min. :-0.0005417
## 1st Qu.:0.9429 1st Qu.:0.7358 1st Qu.: 0.1294366
## Median :0.9577 Median :0.7757 Median : 0.1850655
## Mean :0.9559 Mean :0.7745 Mean : 0.1814488
## 3rd Qu.:0.9722 3rd Qu.:0.8148 3rd Qu.: 0.2275641
## Max. :1.0000 Max. :0.9298 Max. : 0.3817387
## Model coefficients
get(modelac)$betas
## 49 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE -0.03062496
## male 0.07229817
## BMI .
## WHR .
## HEIGHT .
## MASS .
## HG 0.10580374
## PFAT .
## SBP .
## DBP .
## GLY -0.23014855
## HDL .
## LDL -0.35487543
## Ca 0.06963822
## P -0.23750534
## Mg .
## FE .
## SATTRF .
## HGB .
## MCV 0.20337668
## PTH 0.20244539
## UA .
## FOLATE 0.59144460
## uPCR -0.91171129
## AL_adult -0.48157456
## BREAKS .
## log2_FAT -0.02620384
## log2_FFM .
## log2_TC -0.40941964
## log2_TG .
## log2_Se -0.17852882
## log2_Zn -0.48824639
## log2_TIBC .
## log2_FERR -0.46820241
## log2_TRF .
## log2_TRFINDEX 0.12147493
## log2_STRF .
## log2_CTX .
## log2_P1NP 0.37553616
## log2_UI -0.18251390
## log2_UREA -0.50687940
## log2_CREA -0.41883894
## log2_HOLOTC .
## log2_HCY .
## log2_MMA .
## log2_VIT_D 0.42509516
## log2_uCa -0.01533056
## log2_uCCR -0.14743311
## log2_uICR -0.131910123.2.2.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4736842
## Fit and validate model
modelac <- 'pred_complet_adult_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.2626192
## auc 0.8857778
## auc_optimism_corrected 0.6661926
## auc_optimism_corrected_CIL 0.4875733
## auc_optimism_corrected_CIU 0.8360262
## accuracy 0.8315789
## accuracy_optimism_corrected 0.6081130
## accuracy_optimism_corrected_CIL 0.4201613
## accuracy_optimism_corrected_CIU 0.7670814
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.7145 Min. :0.4091 Min. :-0.05124 Min. :0.5851
## 1st Qu.:0.9428 1st Qu.:0.6067 1st Qu.: 0.22871 1st Qu.:0.8737
## Median :0.9696 Median :0.6685 Median : 0.30101 Median :0.9149
## Mean :0.9628 Mean :0.6662 Mean : 0.29657 Mean :0.9075
## 3rd Qu.:0.9903 3rd Qu.:0.7274 3rd Qu.: 0.36386 3rd Qu.:0.9479
## Max. :1.0000 Max. :0.9226 Max. : 0.57402 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.2917 Min. :-0.1104
## 1st Qu.:0.5556 1st Qu.: 0.2307
## Median :0.6154 Median : 0.2991
## Mean :0.6081 Mean : 0.2994
## 3rd Qu.:0.6667 3rd Qu.: 0.3700
## Max. :0.8125 Max. : 0.5715
## Model coefficients
get(modelac)$betas
## 49 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE .
## male .
## BMI .
## WHR .
## HEIGHT .
## MASS .
## HG .
## PFAT .
## SBP .
## DBP 0.01422657
## GLY .
## HDL .
## LDL .
## Ca .
## P .
## Mg .
## FE .
## SATTRF .
## HGB .
## MCV -0.01252481
## PTH .
## UA .
## FOLATE 0.44187862
## uPCR -0.11239698
## AL_adult .
## BREAKS .
## log2_FAT .
## log2_FFM .
## log2_TC -0.25571695
## log2_TG .
## log2_Se -0.20245555
## log2_Zn -0.20792583
## log2_TIBC .
## log2_FERR -0.11934472
## log2_TRF .
## log2_TRFINDEX 0.28357998
## log2_STRF 0.12923032
## log2_CTX .
## log2_P1NP .
## log2_UI -0.31195970
## log2_UREA .
## log2_CREA .
## log2_HOLOTC -0.20347552
## log2_HCY 0.02254704
## log2_MMA 0.25541403
## log2_VIT_D 0.01330014
## log2_uCa .
## log2_uCCR .
## log2_uICR -0.371243643.2.2.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6715328
## Fit and validate model
modelac <- 'pred_complet_adult_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.6000000
## lambda 0.1065992
## auc 0.8246377
## auc_optimism_corrected 0.6567516
## auc_optimism_corrected_CIL 0.4814422
## auc_optimism_corrected_CIU 0.8026877
## accuracy 0.6934307
## accuracy_optimism_corrected 0.6550346
## accuracy_optimism_corrected_CIL 0.5288291
## accuracy_optimism_corrected_CIU 0.7858469
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.2786 Min. :-0.02315 Min. :0.6496
## 1st Qu.:0.9106 1st Qu.:0.6135 1st Qu.: 0.21380 1st Qu.:0.8102
## Median :0.9407 Median :0.6630 Median : 0.27586 Median :0.8540
## Mean :0.9337 Mean :0.6568 Mean : 0.27691 Mean :0.8514
## 3rd Qu.:0.9633 3rd Qu.:0.7098 3rd Qu.: 0.33902 3rd Qu.:0.8972
## Max. :1.0000 Max. :0.8824 Max. : 0.70024 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.4048 Min. :-0.1458
## 1st Qu.:0.6087 1st Qu.: 0.1284
## Median :0.6545 Median : 0.1922
## Mean :0.6550 Mean : 0.1963
## 3rd Qu.:0.7000 3rd Qu.: 0.2608
## Max. :0.8333 Max. : 0.5588
## Model coefficients
get(modelac)$betas
## 49 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE .
## male .
## BMI .
## WHR .
## HEIGHT .
## MASS .
## HG .
## PFAT .
## SBP .
## DBP -0.4100252897
## GLY .
## HDL .
## LDL -0.0298843600
## Ca .
## P .
## Mg .
## FE .
## SATTRF .
## HGB .
## MCV 0.1445031277
## PTH 0.3573358265
## UA .
## FOLATE .
## uPCR -0.2804439274
## AL_adult -0.1792452846
## BREAKS .
## log2_FAT .
## log2_FFM .
## log2_TC .
## log2_TG .
## log2_Se .
## log2_Zn .
## log2_TIBC .
## log2_FERR .
## log2_TRF .
## log2_TRFINDEX .
## log2_STRF .
## log2_CTX .
## log2_P1NP .
## log2_UI .
## log2_UREA -0.2367771870
## log2_CREA .
## log2_HOLOTC 0.0791558791
## log2_HCY -0.0009791532
## log2_MMA -0.1251544859
## log2_VIT_D .
## log2_uCa .
## log2_uCCR .
## log2_uICR .3.2.3 Reduced model
Baseline + clinical characteristics not affected with vitamin supplementation
3.2.3.1 Prepare data
Open code
dat_pred <- dat_adult_imputed %>%
select(GRP, FAM, AGE, male,
BMI:Mg,
HGB:UA,
AL_adult:log2_Zn,
log2_CTX, log2_P1NP,
log2_UREA:log2_CREA
)
dat_features <- dat_pred %>% select(-FAM, -GRP) %>% as.matrix
colnames(dat_features)
## [1] "AGE" "male" "BMI" "WHR" "HEIGHT" "MASS"
## [7] "HG" "PFAT" "SBP" "DBP" "GLY" "HDL"
## [13] "LDL" "Ca" "P" "Mg" "HGB" "MCV"
## [19] "PTH" "UA" "AL_adult" "BREAKS" "log2_FAT" "log2_FFM"
## [25] "log2_TC" "log2_TG" "log2_Se" "log2_Zn" "log2_CTX" "log2_P1NP"
## [31] "log2_UREA" "log2_CREA"3.2.3.2 VN vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VG') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6478873
## Fit and validate model
modelac <- 'pred_complex_adult_VN_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
get(modelac)$plot
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Open code
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 1.00000000
## lambda 0.05634876
## auc 0.87108696
## auc_optimism_corrected 0.78648076
## auc_optimism_corrected_CIL 0.67192971
## auc_optimism_corrected_CIU 0.89560516
## accuracy 0.77464789
## accuracy_optimism_corrected 0.71226811
## accuracy_optimism_corrected_CIL 0.57883333
## accuracy_optimism_corrected_CIU 0.82352941
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.8386 Min. :0.5562 Min. :-0.05388 Min. :0.7376
## 1st Qu.:0.9251 1st Qu.:0.7496 1st Qu.: 0.10153 1st Qu.:0.8389
## Median :0.9403 Median :0.7892 Median : 0.14923 Median :0.8671
## Mean :0.9389 Mean :0.7865 Mean : 0.15240 Mean :0.8661
## 3rd Qu.:0.9557 3rd Qu.:0.8286 3rd Qu.: 0.20162 3rd Qu.:0.8944
## Max. :1.0000 Max. :0.9120 Max. : 0.43843 Max. :0.9930
## accuracy_resamp_validation accuracy_optimism
## Min. :0.4694 Min. :-0.08059
## 1st Qu.:0.6667 1st Qu.: 0.09333
## Median :0.7177 Median : 0.15256
## Mean :0.7123 Mean : 0.15384
## 3rd Qu.:0.7600 3rd Qu.: 0.21353
## Max. :0.8600 Max. : 0.45925
## Model coefficients
get(modelac)$betas
## 32 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE .
## male .
## BMI .
## WHR .
## HEIGHT .
## MASS .
## HG .
## PFAT .
## SBP .
## DBP .
## GLY .
## HDL .
## LDL .
## Ca .
## P .
## Mg .
## HGB .
## MCV .
## PTH .
## UA .
## AL_adult -0.5110666
## BREAKS .
## log2_FAT .
## log2_FFM .
## log2_TC -1.0544345
## log2_TG .
## log2_Se .
## log2_Zn -0.5595906
## log2_CTX .
## log2_P1NP 0.3181682
## log2_UREA -0.9464709
## log2_CREA -0.11490423.2.3.3 VG vs OM
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'VN') %>%
mutate(outcome = if_else(GRP == 'VG', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.4736842
## Fit and validate model
modelac <- 'pred_complex_adult_VG_OM'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.2000000
## lambda 0.6355181
## auc 0.5000000
## auc_optimism_corrected 0.5050453
## auc_optimism_corrected_CIL 0.3253598
## auc_optimism_corrected_CIU 0.6672917
## accuracy 0.5263158
## accuracy_optimism_corrected 0.4761204
## accuracy_optimism_corrected_CIL 0.2872222
## accuracy_optimism_corrected_CIU 0.6461575
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.2125 Min. :-0.1105 Min. :0.5106
## 1st Qu.:0.7939 1st Qu.:0.4512 1st Qu.: 0.2456 1st Qu.:0.6915
## Median :0.8621 Median :0.5064 Median : 0.3411 Median :0.7696
## Mean :0.8357 Mean :0.5050 Mean : 0.3306 Mean :0.7549
## 3rd Qu.:0.9027 3rd Qu.:0.5663 3rd Qu.: 0.4427 3rd Qu.:0.8316
## Max. :1.0000 Max. :0.7593 Max. : 0.7040 Max. :1.0000
## accuracy_resamp_validation accuracy_optimism
## Min. :0.1250 Min. :-0.1480
## 1st Qu.:0.4133 1st Qu.: 0.1940
## Median :0.4848 Median : 0.2881
## Mean :0.4761 Mean : 0.2788
## 3rd Qu.:0.5385 3rd Qu.: 0.3761
## Max. :0.7059 Max. : 0.6260
## Model coefficients
get(modelac)$betas
## 32 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE 0
## male .
## BMI .
## WHR .
## HEIGHT .
## MASS .
## HG .
## PFAT .
## SBP .
## DBP .
## GLY .
## HDL .
## LDL .
## Ca .
## P .
## Mg .
## HGB .
## MCV .
## PTH .
## UA .
## AL_adult .
## BREAKS .
## log2_FAT .
## log2_FFM .
## log2_TC .
## log2_TG .
## log2_Se .
## log2_Zn .
## log2_CTX .
## log2_P1NP .
## log2_UREA .
## log2_CREA .3.2.3.4 VN vs VG
Open code
## Prepare working data
dat_pred_temp <- dat_pred %>%
filter(GRP != 'OM') %>%
mutate(outcome = if_else(GRP == 'VN', 1, 0)) %>%
select(-GRP)
mean(dat_pred_temp$outcome)
## [1] 0.6715328
## Fit and validate model
modelac <- 'pred_complex_adult_VN_VG'
assign(modelac, run(
clustered_glmnet(
data = dat_pred_temp,
outcome = 'outcome',
clust_id = 'FAM',
sample_method = 'atypboot',
N = 500,
alphas = seq(0, 1, by = 0.2),
seed = 368),
path = paste0('gitignore/run/', modelac),
reuse = TRUE))
## Calibration curve
#get(modelac)$plot
## Model summary
t(get(modelac)$model_summary)
## [,1]
## alpha 0.0000000
## lambda 0.7353789
## auc 0.8140097
## auc_optimism_corrected 0.6556023
## auc_optimism_corrected_CIL 0.4932692
## auc_optimism_corrected_CIU 0.7907299
## accuracy 0.7080292
## accuracy_optimism_corrected 0.6500979
## accuracy_optimism_corrected_CIL 0.5185185
## accuracy_optimism_corrected_CIU 0.7659574
get(modelac)$valid_performances %>% summary
## auc_resamp_test auc_resamp_validation auc_optimism accuracy_resamp_test
## Min. :0.5000 Min. :0.3782 Min. :-0.03536 Min. :0.6475
## 1st Qu.:0.8682 1st Qu.:0.6118 1st Qu.: 0.17729 1st Qu.:0.7704
## Median :0.9001 Median :0.6599 Median : 0.23573 Median :0.8058
## Mean :0.8931 Mean :0.6556 Mean : 0.23749 Mean :0.8074
## 3rd Qu.:0.9261 3rd Qu.:0.7059 3rd Qu.: 0.29642 3rd Qu.:0.8456
## Max. :0.9897 Max. :0.8633 Max. : 0.58373 Max. :0.9485
## accuracy_resamp_validation accuracy_optimism
## Min. :0.4048 Min. :-0.08852
## 1st Qu.:0.6042 1st Qu.: 0.09059
## Median :0.6522 Median : 0.14748
## Mean :0.6501 Mean : 0.15732
## 3rd Qu.:0.6949 3rd Qu.: 0.21683
## Max. :0.8077 Max. : 0.46158
## Model coefficients
get(modelac)$betas
## 32 x 1 sparse Matrix of class "dgCMatrix"
## s0
## AGE -0.130969640
## male 0.046847677
## BMI -0.025892148
## WHR -0.071845276
## HEIGHT -0.057146386
## MASS -0.052079055
## HG 0.018645667
## PFAT -0.072209270
## SBP 0.025161999
## DBP -0.218387885
## GLY -0.054679974
## HDL 0.007689509
## LDL -0.110257424
## Ca 0.016184827
## P -0.086481596
## Mg 0.011665361
## HGB 0.007370975
## MCV 0.179203214
## PTH 0.210230420
## UA 0.057920958
## AL_adult -0.204033621
## BREAKS -0.084629183
## log2_FAT -0.065394397
## log2_FFM -0.007507923
## log2_TC -0.078258475
## log2_TG 0.017640442
## log2_Se 0.003677581
## log2_Zn -0.139279629
## log2_CTX 0.105505924
## log2_P1NP 0.094490749
## log2_UREA -0.199522992
## log2_CREA -0.0781046813.3 Results summaries
3.3.1 Prepare data
Open code
models_names <- c(
"pred_baseline_child_all_VN_OM",
"pred_baseline_child_all_VG_OM",
"pred_baseline_child_all_VN_VG",
"pred_complex_child_all_VN_OM",
"pred_complex_child_all_VG_OM",
"pred_complex_child_all_VN_VG",
"pred_complet_child_all_VN_OM",
"pred_complet_child_all_VG_OM",
"pred_complet_child_all_VN_VG",
"pred_baseline_adult_VN_OM",
"pred_baseline_adult_VG_OM",
"pred_baseline_adult_VN_VG",
"pred_complex_adult_VN_OM",
"pred_complex_adult_VG_OM",
"pred_complex_adult_VN_VG",
"pred_complet_adult_VN_OM",
"pred_complet_adult_VG_OM",
"pred_complet_adult_VN_VG"
)
tr <- get(models_names[1])$model_summary[
c(
"auc_optimism_corrected",
"auc_optimism_corrected_CIL",
"auc_optimism_corrected_CIU"
)
]
for (i in 2:length(models_names)) {
tr <- dplyr::bind_rows(tr, get(models_names[i])$model_summary[
c(
"auc_optimism_corrected",
"auc_optimism_corrected_CIL",
"auc_optimism_corrected_CIU"
)
])
}
tr <- data.frame(models_names,
y_pos = c(
rep(0.9, 3),
rep(0.8, 3),
rep(0.7, 3),
rep(0.5, 3),
rep(0.4, 3),
rep(0.3, 3)
),
age_group = c(
rep("children", 9),
rep("adults", 9)
),
contrast = str_sub(models_names, -5, -1),
model = str_sub(models_names, 6, 12),
tr
) %>%
mutate(
model = fct_recode(model,
"baseline" = "baselin",
"reduced" = "complex",
"full" = "complet"
)
) %>%
rename(
AUC = auc_optimism_corrected,
CI_L = auc_optimism_corrected_CIL,
CI_U = auc_optimism_corrected_CIU
) %>%
select(
age_group, model, contrast, AUC, CI_L, CI_U, models_names, y_pos
)3.3.2 Models performance summary
Open code
tr %>%
mutate(`AUC [95% CI]` = paste0(
round(AUC, 2),
" [",
round(CI_L, 2),
", ",
round(CI_U, 2),
"]"
)) %>%
select(age_group:contrast, `AUC [95% CI]`) %>%
kableExtra::kable(caption = "Results of elastic net predictive models distinguishing between diet groups. `age_group` indicates the cohort (children or adults). `model` specifies the covariates included: `baseline` includes only basic patient characteristics (age, sex, breastfeeding-related covariates and birth weight for children); `reduced` adds clinical characteristics not primarily driven by supplementation; `full` includes all clinical characteristics. `contrast` denotes the pair of diet groups used for model training: `VN_OM` (vegans vs. omnivores), `VG_VN` (vegetarians vs. vegans), and `OM_VG` (omnivores vs. vegetarians). `AUC [95% CI]` represents the out-of-sample AUC estimated via clustered bootstrap, with corresponding 95% confidence intervals.")| age_group | model | contrast | AUC [95% CI] |
|---|---|---|---|
| children | baseline | VN_OM | 0.54 [0.36, 0.69] |
| children | baseline | VG_OM | 0.57 [0.33, 0.78] |
| children | baseline | VN_VG | 0.59 [0.39, 0.76] |
| children | reduced | VN_OM | 0.65 [0.46, 0.8] |
| children | reduced | VG_OM | 0.42 [0.21, 0.64] |
| children | reduced | VN_VG | 0.64 [0.47, 0.79] |
| children | full | VN_OM | 0.86 [0.74, 0.95] |
| children | full | VG_OM | 0.59 [0.32, 0.85] |
| children | full | VN_VG | 0.7 [0.52, 0.88] |
| adults | baseline | VN_OM | 0.54 [0.38, 0.71] |
| adults | baseline | VG_OM | 0.47 [0.3, 0.59] |
| adults | baseline | VN_VG | 0.55 [0.41, 0.73] |
| adults | reduced | VN_OM | 0.79 [0.67, 0.9] |
| adults | reduced | VG_OM | 0.51 [0.33, 0.67] |
| adults | reduced | VN_VG | 0.66 [0.49, 0.79] |
| adults | full | VN_OM | 0.85 [0.75, 0.94] |
| adults | full | VG_OM | 0.67 [0.49, 0.84] |
| adults | full | VN_VG | 0.66 [0.48, 0.8] |
3.3.3 Comparison of baseline vs more complex models
Open code
baseline <- tr %>% filter(model == "baseline")
full <- tr %>% filter(model == "full")
reduced <- tr %>% filter(model == "reduced")
compmod <- compare_predmod(
get(baseline$models_names[1]),
get(reduced$models_names[1])
)
comp_table <- data.frame(
age_group = baseline$age_group[1],
contrast = baseline$contrast[1],
auc_gain = compmod[1],
CI_L = compmod[2],
CI_U = compmod[3],
P = compmod[4],
models = interaction(baseline$model[1], reduced$model[1])
)
for (i in 2:nrow(baseline)) {
compmod <- compare_predmod(
get(baseline$models_names[i]),
get(reduced$models_names[i])
)
comp_table_tmp <- data.frame(
age_group = baseline$age_group[i],
contrast = baseline$contrast[i],
auc_gain = compmod[1],
CI_L = compmod[2],
CI_U = compmod[3],
P = compmod[4],
models = interaction(baseline$model[i], reduced$model[i])
)
comp_table <- dplyr::bind_rows(comp_table, comp_table_tmp)
i <- i + 1
}
comp_table
## age_group contrast auc_gain CI_L CI_U P
## 1 children VN_OM 0.11529848 -0.09366843 0.2846901 0.23752495
## 2 children VG_OM -0.15032444 -0.41028705 0.1104278 0.24151697
## 3 children VN_VG 0.05133982 -0.15804438 0.2717500 0.55289421
## 4 adults VN_OM 0.24414073 0.04718529 0.4171286 0.01796407
## 5 adults VG_OM 0.03800912 -0.17700521 0.2350962 0.68063872
## 6 adults VN_VG 0.10676236 -0.12723077 0.2916598 0.32135729
## models
## 1 baseline.reduced
## 2 baseline.reduced
## 3 baseline.reduced
## 4 baseline.reduced
## 5 baseline.reduced
## 6 baseline.reduced
dim1 <- nrow(comp_table)
for (i in (dim1 + 1):(dim1 + nrow(baseline))) {
compmod <- compare_predmod(
get(baseline$models_names[i - dim1]),
get(full$models_names[i - dim1])
)
comp_table_tmp <- data.frame(
age_group = baseline$age_group[i - dim1],
contrast = baseline$contrast[i - dim1],
auc_gain = compmod[1],
CI_L = compmod[2],
CI_U = compmod[3],
P = compmod[4],
models = interaction(baseline$model[i - dim1], full$model[i - dim1])
)
comp_table <- dplyr::bind_rows(comp_table, comp_table_tmp)
}
comp_table %>%
mutate(
`AUC gain [95% CI]` = paste0(
round(auc_gain, 2), " [", round(CI_L, 2), ", ", round(CI_U, 2), "]"
),
models = fct_recode(models,
"reduced vs baseline" = "baseline.reduced",
"full vs baseline" = "baseline.full"
)
) %>%
select(
models, age_group, contrast, `AUC gain [95% CI]`, P
) %>%
kableExtra::kable(
caption = "Comparison of elastic net predictive models evaluating the gain in discriminative performance when including additional clinical characteristics. The `models` column compares either `reduced vs baseline` (adding clinical characteristics not primarily driven by supplementation) or `full vs baseline` (adding all clinical characteristics). `age_group` indicates the cohort (children or adults). `contrast` specifies the diet groups used for model training: `VN_OM` (vegans vs. omnivores), `VG_OM` (vegetarians vs. omnivores), and `VN_VG` (vegans vs. vegetarians). `AUC gain [95% CI]` represents the estimated improvement in out-of-sample AUC with a more complex model, with clustered bootstrap providing corresponding 95% confidence intervals. `P` indicates the p-value assessing the statistical significance of the AUC gain."
)| models | age_group | contrast | AUC gain [95% CI] | P |
|---|---|---|---|---|
| reduced vs baseline | children | VN_OM | 0.12 [-0.09, 0.28] | 0.2375250 |
| reduced vs baseline | children | VG_OM | -0.15 [-0.41, 0.11] | 0.2415170 |
| reduced vs baseline | children | VN_VG | 0.05 [-0.16, 0.27] | 0.5528942 |
| reduced vs baseline | adults | VN_OM | 0.24 [0.05, 0.42] | 0.0179641 |
| reduced vs baseline | adults | VG_OM | 0.04 [-0.18, 0.24] | 0.6806387 |
| reduced vs baseline | adults | VN_VG | 0.11 [-0.13, 0.29] | 0.3213573 |
| full vs baseline | children | VN_OM | 0.33 [0.16, 0.51] | 0.0019960 |
| full vs baseline | children | VG_OM | 0.02 [-0.28, 0.35] | 0.9281437 |
| full vs baseline | children | VN_VG | 0.11 [-0.12, 0.37] | 0.3373253 |
| full vs baseline | adults | VN_OM | 0.31 [0.12, 0.51] | 0.0019960 |
| full vs baseline | adults | VG_OM | 0.2 [-0.01, 0.42] | 0.0618762 |
| full vs baseline | adults | VN_VG | 0.11 [-0.12, 0.31] | 0.3173653 |
3.3.4 Visualize
Open code
p1 <- tr %>%
mutate(
contrast = factor(contrast, levels = c("VN_OM", "VG_OM", "VN_VG")),
model = factor(model, levels = c("baseline", "reduced", "full"))
) %>%
ggplot(aes(x = AUC, y = model, color = model)) +
geom_point(size = 5, shape = 18) +
geom_segment(
aes(
x = CI_L, xend = CI_U,
y = model, yend = model
),
linewidth = 1.1
) +
coord_cartesian(xlim = c(0.1, 1)) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
facet_wrap(~ age_group + contrast,
labeller = labeller(
contrast = c(
"VN_OM" = "VN vs OM",
"VG_OM" = "VG vs OM",
"VN_VG" = "VN vs VG"
)
)
) +
theme(
axis.title.y = element_blank(),
legend.position = "bottom"
) +
labs(
x = "Out-of-sample AUC",
color = "Model"
) +
scale_color_manual(values = c(
"baseline" = "grey50",
"full" = "red",
"reduced" = "purple"
))
p1
VN vs OM, left), vegetarians vs. omnivores (VG vs OM, middle), and vegans vs. vegetarians (VN vs VG, right). Models include baseline (only basic patient characteristics: age, sex, breastfeeding-related covariates and birth weight for children), reduced (adding clinical characteristics not primarily driven by supplementation), and full (including all clinical characteristics). Points represent out-of-sample AUC estimates, indicating the ability of models to discriminate between diet groups, with lines showing 95% confidence intervals estimated via clustered bootstrap.4 Reproducibility
Open code
sessionInfo()
## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=cs_CZ.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=cs_CZ.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=cs_CZ.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=cs_CZ.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Prague
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mice_3.17.0 patchwork_1.2.0 ggrepel_0.9.5 robustlmm_3.3-1
## [5] gridExtra_2.3 pheatmap_1.0.12 performance_0.12.2 quantreg_5.98
## [9] SparseM_1.81 bayesplot_1.8.1 ggdist_3.3.2 kableExtra_1.4.0
## [13] lubridate_1.8.0 corrplot_0.92 arm_1.12-2 MASS_7.3-64
## [17] projpred_2.0.2 glmnet_4.1-8 boot_1.3-31 cowplot_1.1.1
## [21] pROC_1.18.0 mgcv_1.9-1 nlme_3.1-167 openxlsx_4.2.5
## [25] flextable_0.9.6 sjPlot_2.8.16 car_3.1-2 carData_3.0-5
## [29] gtsummary_2.0.2 emmeans_1.10.4 ggpubr_0.4.0 lme4_1.1-35.5
## [33] Matrix_1.7-0 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [37] purrr_1.0.2 readr_2.1.2 tidyr_1.3.1 tibble_3.2.1
## [41] ggplot2_3.5.1 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.4.3 later_1.3.0 gamm4_0.2-6
## [4] cellranger_1.1.0 datawizard_0.12.2 rpart_4.1.24
## [7] reprex_2.0.1 lifecycle_1.0.4 rstatix_0.7.0
## [10] lattice_0.22-5 insight_0.20.2 backports_1.5.0
## [13] magrittr_2.0.3 rmarkdown_2.27 yaml_2.3.5
## [16] httpuv_1.6.5 skimr_2.1.5 zip_2.2.0
## [19] askpass_1.1 DBI_1.1.2 minqa_1.2.4
## [22] RColorBrewer_1.1-2 multcomp_1.4-18 abind_1.4-5
## [25] rvest_1.0.2 nnet_7.3-20 TH.data_1.1-0
## [28] sandwich_3.0-1 gdtools_0.3.7 crul_1.5.0
## [31] MatrixModels_0.5-3 svglite_2.1.3 codetools_0.2-19
## [34] xml2_1.3.3 tidyselect_1.2.1 shape_1.4.6
## [37] farver_2.1.0 ggeffects_1.7.0 httpcode_0.3.0
## [40] base64enc_0.1-3 matrixStats_1.3.0 jsonlite_1.8.8
## [43] mitml_0.4-3 ellipsis_0.3.2 ggridges_0.5.3
## [46] survival_3.7-0 iterators_1.0.14 systemfonts_1.0.4
## [49] foreach_1.5.2 tools_4.4.3 ragg_1.2.1
## [52] Rcpp_1.0.13 glue_1.7.0 pan_1.6
## [55] xfun_0.46 distributional_0.4.0 loo_2.4.1
## [58] withr_3.0.1 fastmap_1.2.0 fansi_1.0.6
## [61] openssl_1.4.6 digest_0.6.37 R6_2.5.1
## [64] mime_0.12 estimability_1.5.1 textshaping_0.3.6
## [67] colorspace_2.0-2 utf8_1.2.4 generics_0.1.3
## [70] fontLiberation_0.1.0 data.table_1.15.4 robustbase_0.93-9
## [73] httr_1.4.2 htmlwidgets_1.6.4 pkgconfig_2.0.3
## [76] gtable_0.3.0 htmltools_0.5.8.1 fontBitstreamVera_0.1.1
## [79] scales_1.3.0 knitr_1.48 rstudioapi_0.16.0
## [82] tzdb_0.2.0 uuid_1.0-3 coda_0.19-4
## [85] curl_4.3.2 nloptr_2.0.0 repr_1.1.7
## [88] zoo_1.8-9 sjlabelled_1.2.0 parallel_4.4.3
## [91] pillar_1.9.0 vctrs_0.6.5 promises_1.2.0.1
## [94] jomo_2.7-3 dbplyr_2.1.1 xtable_1.8-4
## [97] evaluate_1.0.0 fastGHQuad_1.0.1 mvtnorm_1.1-3
## [100] cli_3.6.3 compiler_4.4.3 rlang_1.1.4
## [103] crayon_1.5.0 rstantools_2.1.1 ggsignif_0.6.3
## [106] labeling_0.4.2 modelr_0.1.8 plyr_1.8.6
## [109] sjmisc_2.8.10 fs_1.6.4 stringi_1.7.6
## [112] viridisLite_0.4.0 assertthat_0.2.1 munsell_0.5.0
## [115] fontquiver_0.2.1 sjstats_0.19.0 hms_1.1.1
## [118] gfonts_0.2.0 shiny_1.9.1 haven_2.4.3
## [121] broom_1.0.6 DEoptimR_1.0-10 readxl_1.3.1
## [124] officer_0.6.6